This work proposes a novel smoothing technique for functional data observed over space and time, that penalizes the misfit with respect to a timedependent nonlinear partial differential equation. This method extends existing physics-informed statistical models, so far explored in the case of penalty terms involving linear partial differential equations. Their use in the regularizing term of a smoothing model broaden extensively the scope of this method, but introduces new challenges in both theory and numerical implementation. Through a simulation study, we show that the proposed method outperforms simpler models, which are unable to incorporate complex physics information in the statistical model. These results demonstrate the potential of combining advanced physical knowledge with statistical modelling, paving the way for future applications to real-world data, as briefly discussed in the context of neuroimaging data.