I would like to solve 3 differential equations for 3 unknowns. So I wrote a MATLAB code, which solves (using the '\' operator) these equation using a linear system of equations (in which the 3 equations are coupled). So the linear system has a size of (ny*nx*3,ny*nx*3) [nx/ny are the # of nodes in x/y - direction]. I used a finite-difference method.
Solved for $F$:
$\nabla u - \nabla a \nabla F = \nabla a B_1$
where $u$ is the unknown of a different equation, and $a$, $\alpha$, $t$ and $B$ are parameters.
The parameter $a$ produces very different results of $F$, as shown in the figure below [dotted-black line-$a$ parameter, which is here normalized to the range of $F$; solid red line - solution of $F$]. These results seem to be correct, since they reflect the solution of a benchmark.
Now I would like to isolate these equations and pass the solution of one equation to the the other one (as a parameter on the right-hand-side of the equation). So I store the solution of $u$ from the previous ('coupled') approach and pass it as a parameter to the equation, which is written above. Again I use MATLAB, finite differences and the '\' operator. But this time the linear system of equation only has a size of (ny * nx,ny * nx), since I solve only for one unknown ($F$). I expect the solution of this isolated system to be very similar to the one of the coupled system, since I solve the same equation and have the same values for $u$, $a$, $B$, spatial step size ... !
Unfortunately, the results are very different from the previous ('coupled') approach, if $a$ is very low (see the solid blue line on the lower figure). Especially where $\nabla a$ becomes very low, there is a great discrepancy. Why is such an approach numerical instable? Is it because of floating point inaccuracy? How could I make the problem work?
