To solve the split between structural recovery and continuous dynamics’ characterisation in gene regulatory network identification, we put forward a Neural ODE-GRN joint modelling method which is based on neural ordinary differential equations. This method uses time-series bulk RNA-seq and pseudo-time-series single-cell RNA-seq expression matrices as input materials, introduces continuous-time evolution equations into a hidden state space, and combines sparse structural masks, symbolic constraints, biological priors and temporal consistency loss to realize co-optimisation of regulatory adjacency matrix learning and expression trajectory reconstruction. Therefore, validation works were conducted with simulated data from GeneNetWeaver and SERGIO, as well as real-world data from hESCs and mESCs. Hence, results show that on a 100-gene simulated network, the model’s AUPR achieved 0.631, which is an improvement of 0.042 to 0.114 compared with various baselines; under 20% noise conditions, the AUROC kept above 0.80; on real-world data, the trajectory correlation coefficient reached 0.914-0.928, at the same time maintaining good topological consistency for hub nodes and representative regulatory modules. Therefore, this method provides a new computational tool for explaining dynamic regulatory mechanisms in developmental differentiation and disease progression.