The advanced single-cell RNA sequencing (scRNA-seq) technology allows us to measure the temporal dynamics of gene expression over multiple time points to gain an understanding of previously unknown biological diversity. However, there is currently a lack of efficient computational tools tailored for scRNA-seq data analysis, which can simultaneously analyze the entire sequence across different time points as well as accounting for temporal batch effects. In this paper, we present a non-parametric statistical method called TDEseq that takes full advantage of smoothing splines basis functions to account for the dependence of multiple time points, and uses hierarchical structure linear additive mixed models to model the correlated cells within an individual. As a result, TDEseq demonstrates powerful performance in identifying four potential temporal expression patterns within a specific cell type, including growth, recession, peak, and trough. Extensive simulation studies and the analysis of four published scRNA-seq datasets show that TDEseq can produce well-calibrated p-values and up to 20% power gain over the existing methods for detecting temporal gene expression patterns, even with the case in large heterogeneous. TDEseq is also capable of handling unwanted confounding factors that may be hidden in biological processes, thereby enabling advancements in investigations that utilize time-resolved or time-course scRNA-seq data, particularly in the multi-sample multi-stage studies.