Calculating the jacobian matrix for a time-dependent system

Question

For testing purposes I would like to calculate the Jacobian matrix of the time-dependent heat equation, i.e. $$\partial_tu=\nabla^2u+f$$ As far as I understood, the Jacobian can be calculated using $$F(u)=\partial_tu-\nabla^2u-f$$ and (with $\delta u$ a small step) $$\begin{split}F'(u)(\delta u)&=\lim\limits_{h\rightarrow0}\frac{F(u+h\delta u)-F(u)}{h}\\ &=\lim\limits_{h\rightarrow0}\frac{\partial_t(u+h\delta u)-\nabla^2(u+h\delta u)-f-\partial_tu+\nabla^2u+f}{h}\\ &=\lim\limits_{h\rightarrow0}\frac{h\partial_t(\delta u)-h\nabla^2(\delta u)}{h}\\ &=\partial_t(\delta u)-\nabla^2(\delta u)\forall h\neq0\end{split}$$ That result can be rewritten into $$F'(u)(\delta u)=\frac{\delta u_1-\delta u_0}{dt}-\nabla^2(\delta u)$$ Is that correct? If yes, what would $\delta u_1$ and $\delta u_0$ be? Or is there a mistake in the last step?