Finite differences and Neumann boundary conditions
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 & & ... 935