I am dealing with a highly nonlinear system of two PDEs. I already have a code to solve the system in case of Dirichlet boundary conditions. The explicit system is:
$$ \begin{eqnarray*} \partial_{t}u & = & D_{11}\partial_{x}^{2}u+\partial_{x}D_{11}\cdot\partial_{x}u+D_{12}\partial_{x}^{2}v+\partial_{x}D_{12}\cdot\partial_{x}v\\ \partial_{t}v & = & D_{21}\partial_{x}^{2}u+\partial_{x}D_{21}\cdot\partial_{x}u+D_{22}\partial_{x}^{2}v+\partial_{x}D_{22}\cdot\partial_{x}v \end{eqnarray*}$$
with boundary conditions:
$$u(0,t)=0 \ (\partial_x v)(0,t) = 0$$
And the discretized form of the scheme are:
$$\begin{eqnarray*} -D_{21}(x_{j})u_{n+1,j-1}-D_{22}(x_{j})v_{n+1,j-1}\\ {}[\mu+2D_{11}(x_{j})]u_{n+1,j}+2D_{12}(x_{j})v_{n+1,j}\\ -D_{11}(x_{j})u_{n+1,j+1}-D_{12}(x_{j})v_{n+1,j+1} & = & \frac{1}{4}[4\mu-D_{11}(x_{j+1})+D_{11}(x_{j-1})]u_{n,j}\\ & & -[D_{12}(x_{j+1})-D_{12}(x_{j})]v_{n,j}\\ & & +\frac{1}{4}[D_{11}(x_{j+1})-D_{11}(x_{j-1})]u_{n,j+1}\\ & & +\frac{1}{4}[D_{12}(x_{j+1})-D_{12}(x_{j-1})]v_{n,j+1} \end{eqnarray*}$$
$$ \begin{eqnarray*} -D_{21}(x_{j})u_{n+1,j-1}-D_{22}(x_{j})v_{n+1,j-1}\\ 2D_{21}(x_{j})u_{n+1,j}+[\mu+2D_{22}(x_{j})]v_{n+1,j}\\ -D_{21}(x_{j})u_{n+1,j+1}-D_{22}(x_{j})v_{n+1,j+1} & = & \frac{1}{4}[4\mu-D_{21}(x_{j+1})+D_{21}(x_{j-1})]u_{n,j}\\ & & -[D_{22}(x_{j+1})-D_{22}(x_{j})]v_{n,j}\\ & & +\frac{1}{4}[D_{21}(x_{j+1})-D_{21}(x_{j-1})]u_{n,j+1}\\ & & +\frac{1}{4}[D_{22}(x_{j+1})-D_{22}(x_{j-1})]v_{n,j+1} \end{eqnarray*} $$
That is, I am solving a linear system of the form:
$$A\boldsymbol{x}=\boldsymbol{b}$$
where $A$ has $2J\times 2J$ entries, according to the structure:
$$\left(\begin{array}{cccccccccc} \mbox{1st equation}\rightarrow & | & c_{u}^{j} & c_{v}^{j} & c_{u}^{j+1} & c_{v}^{j+1} & 0 & & ... & 0\\ \mbox{2nd equation}\rightarrow & | & d_{u}^{j} & d_{v}^{j} & d_{u}^{j+1} & d_{v}^{j+1} & 0 & & ... & 0\\ \vdots & | & . & . & . & . & . & . & . & .\\ \mbox{1st}\rightarrow & | & 0 & 0 & c_{u}^{j-1} & c_{v}^{j-1} & c_{u}^{j} & c_{v}^{j} & c_{u}^{j+1} & c_{v}^{j+1}\\ \mbox{2nd }\rightarrow & | & 0 & 0 & d_{u}^{j-1} & d_{v}^{j-1} & d_{u}^{j} & d_{v}^{j} & d_{u}^{j+1} & d_{v}^{j+1}\vdots \end{array}\right)$$
where the $c$'s and $d$'s are the the coefficients corresponding to the scheme above.
Now, I am trying to implement the Neumann boundary condition (the one on $v$), using the ghost point method as discussed, for instance, here for a simpler equation.
Thus, I would like to change my code in order to add some extra terms to the array, but I can't work out what to put in them. Also, in the respective components of the matrix $A$. Because when I write the condition:
$$\frac{u_1-u_{-1}}{(\Delta x)^2}=0$$ (and same for $v$)
I need to eliminate $u_{-1}$ and $v_{-1}$, using also the third and fourth equations I wrote in this post. If I inverted for $u_{-1}$ and $v_{-1}$ (assuming I know how to), I do not see how I could write down a code which works for any kind of scheme I have, without putting the explicit inversion in the matrix.
