# Should the Jacobian of a system of PDEs be calculated from the main equations of the discretised equation?

I am solving a coupled system of non-linear PDEs in 1D. Something like,

$$u_t = F_1(u,v,w) \\ v_t = F_2(u,v,w) \\ w_t = F_3(u,v,w)$$

where each variable is a function of $x$ (the spatial dimension).

So for each point in space I wish to calculate the Jacobian of the vector,

$$\boldsymbol{F} = \pmatrix{F_1(u,v,w) \\ F_2(u,v,w) \\ F_3(u,v,w)}$$

Question

Should I use the PDEs themselves as a starting point to calculate the Jacobian, or is it better to use the discretised form of the equations?

I ask because the equations are very unstable without using exponential fitting, so the discretised form is quite non-trivial. My intuition is that Jacobian should follow from the discretised form as this seems more consistent. Is there a standard approach?

For example, if I were to start with the discretised form, I would simply partially differentiate each term with respect to the appropriate variable.