Say we have a discretised a coupled nonlinear system of two PDEs to give a system of ODEs which approximates the original system, $$ \frac{\partial u}{\partial t} = F_1(t,\boldsymbol{u,v}) \\ \frac{\partial v}{\partial t} = F_2(t,\boldsymbol{u,v}) \\ $$ where $\boldsymbol{u} = [u_1 \cdots u_N]^{T}$ and $\boldsymbol{v} = [v_1 \cdots v_N]^{T}$ are a vector of solution variables with the same number of elements as discretised points in the system. $F$ is a vector which depends on how the PDEs has been discretised.
When PDEs are left in this semi-discrete form (i.e. the transient term is not discretised) they can be solved using conventional ODE solvers, this technique is known as the Method of Lines (MOL).
Normally the ODE solver can be used to step the system forward in time until a desired time point or until the system reaches steady-state. However, can the MOL technique be used to solve directly for the steady-state value?
Internally the ODE solver will be computing a Jacobian and using Newton's method to solve the system, these are the prerequisites for solving the system directly for steady state. So it would seem that it has all the necessary information to do so. Moreover, to solve the above system for the steady state one would apply Newton's method but would set the transient term to zero,
$$ 0 = F_1(t,\boldsymbol{u,v}) \\ 0 = F_2(t,\boldsymbol{u,v}) \\ $$
Is there a way to do this with MOL solvers (e.g. in MATLAB or Python)?