I have a system of 2 non-linear, coupled PDEs that I would like to transform to a stiff ODE system to solve them using the method of lines. The 2 equations:
$$\begin{align} \frac{\partial \phi}{\partial t} &= - \phi^m \frac{P}{R - H(P)(R-1)}\\ \phi^m \frac{P}{R - H(P)(R-1)} &= \nabla \cdot [\phi^n(\nabla P + u_z)] \end{align} $$
$\phi$ and $P$ are my two unknowns, $R$ is a constant, $H(P)$ is the heaviside function of $P$ and $u_z$ is an upward directed unit vector.
The problem is that my second equation is an elliptic equation, and hence I have no time derivative in there. My idea is to transform them to Differential Algebraic Equations (DAE) using finite differences for the space derivatives and then to do index reduction to obtain an ODE system.
I have no experience with DAE, so I am stuck there. If we work in 1D in the z direction, I obtain with finite differences:
$$ \frac{\partial \phi_i}{\partial t} = - \phi_i^m \frac{P_i}{R - H(P_i)(R-1)} \\ (\frac{\phi_i^m \Delta z}{R - H(P_i)(R-1)} + \frac{\phi^n_{i+\frac{1}{2}}}{\Delta z} + \frac{\phi^n_{i-\frac{1}{2}}}{\Delta z}) P_i - \frac{\phi^n_{i-\frac{1}{2}}}{\Delta z} P_{i-1} - \frac{\phi^n_{i+\frac{1}{2}}}{\Delta z} P_{i+1}= (\phi_{i+\frac{1}{2}}^n - \phi_{i-\frac{1}{2}}^n) $$
From what I read, it seems similar to a Hessenberg index-1:
$$\begin{align} y' &= f(t, y, z) \\ 0 &= g(t, y, z) \end{align} $$
This would only require 1 additional derivative.
May I have a bit of help on how to start from there? Pointers to do that using symbolic computation like Sympy could also work!