Nonlinear coupled BVP using relaxation method (finite-difference)

Question

It's my first time posting here and it seems like latex is not working, so I posted some part as an image.

Reference for Relaxation method: Computational Physics by Veseley

---

I have employed the two discussions above in my given nonlinear coupled ODE. The key parameters of the code are the number of grid points n, number of iteration m, and the relaxation parameter w.

It's easier to obtain a residual tolerance that is of order O(10^{-15}) when considering only Dirichlet conditions on both sides. The problems I encountered are as follows,

• There are a few choices of the parameters n, m, and w that lead to a result that does not produce an error, i.e. badly conditioned matrix. However, if I changed a little bit the constants of my setup, say u(a), or even the parameters n and m themselves to get more accuracy, then I need to find another, say w for fixed n and m that leads to a result without error. In general, I could say that my code is not stable given the behavior.
• Given that the code works, the Neumann boundary condition is causing an issue since I cannot get a residual tolerance that is small enough like if the Neumann boundary condition were a Dirichlet condition such that it satisfies both the governing equation and the boundary condition.

I have written my code in Mathematica and displayed a working choice.

(*Equation setup*)
Needs["VariationalMethods`"]
u0 = 10;
m = (1/u[x] - 1/u0);
f = 1 - m z[x]^(d + 1);
L = (Sqrt[-f u'[x]^2 - 2 u'[x] z'[x] + 1]/z[x]^d);(*Lagrangian*)
eq1 = EulerEquations[L, u[x], x];(*Euler-Lagrange equations*)
eq2 = EulerEquations[L, z[x], x];
s = Solve[{eq1, eq2} /. d -> 3, {u''[x], z''[x]}][[1]] // Simplify;
eq01 = (u''[x] - s[[1, 2]]);(*eq01 & eq02 are the nonlinear coupled ODEs*)
eq02 = (z''[x] - s[[2, 2]]);

(*Finite difference and construction of the residual*)
n = 800;(*number of grid points*)
h = (b - a)/(n - 1);(*Step size*)
a = 10^-1;(*a & b are the domain*)
b = 5;
ui = 5;(*ui, zi, are the Dirichlet boundary conditions*)
zi = 10^-1;
uf = 9;(*uf, zf, are the test values for the final grid points u[n], z[n]*)
zf = 5;
up = 0;(*up & zp are the Neumann boundary conditions which is incorporated into the residual down below*)
zp = 0;
rule = {u''[x] -> ((u[i + 1] - 2 u[i] + u[i - 1])/h^2), u'[x] -> ((u[i + 1] - u[i - 1])/(2 h)), u[x] -> u[i], z''[x] -> ((z[i + 1] - 2 z[i] + z[i - 1])/h^2), z'[x] -> ((z[i + 1] - z[i - 1])/(2 h)), z[x] -> z[i]};(*finite difference rule for the interior points*)
resid1 = h^2 eq01 /. rule;(*resid1 & resid2 are the interior points' residuals*)
resid2 = h^2 eq02 /. rule;
residbound1 = (-u[n - 2] + 8 u[n - 1] - 7 u[n] + 6 up h) - 2 h^2 s[[1, 2]] /. {d -> 3, u[x] -> u[n], z[x] -> z[n], u'[x] -> up, z'[x] -> zp};(*residbound1 & residbound2 are the residuals at grid point n*)
residbound2 = (-z[n - 2] + 8 z[n - 1] - 7 z[n] + 6 zp h) - 2 h^2 s[[2, 2]] /. {d -> 3, u[x] -> u[n], z[x] -> z[n], u'[x] -> up, z'[x] -> zp};

(*Construction of sparse matrix*)
parresid1 = {D[resid1, u[i - 1]], D[resid1, z[i - 1]], D[resid1, u[i]], D[resid1, z[i]], D[resid1, u[i + 1]], D[resid1, z[i + 1]]};(*partial derivative of the residual*)
parresid2 = {D[resid2, u[i - 1]], D[resid2, z[i - 1]], D[resid2, u[i]], D[resid2, z[i]], D[resid2, u[i + 1]], D[resid2, z[i + 1]]};
parresidbound1 = D[residbound1, {{u[n - 2], z[n - 2], u[n - 1], z[n - 1], u[n], z[n]}}];(*partial derivative of the boundary residual*)
parresidbound2 = D[residbound2, {{u[n - 2], z[n - 2], u[n - 1], z[n - 1], u[n], z[n]}}];
mat = {parresid1, parresid2};
sparseresidual = Normal[SparseArray[Table[Band[{2 (i - 2) + 1, 2 (i - 2) + 1}] -> {mat}, {i, 2, n - 1}]]];
sparse = Join[{Join[{1}, ConstantArray[0, 2 n - 1]]}, {Join[{0, 1}, ConstantArray[0, 2 n - 2]]}, sparseresidual, {Join[ConstantArray[0, 2 n - 6], parresidbound1]}, {Join[ConstantArray[0, 2 n - 6], parresidbound2]}];(*sparse matrix*)

(*Solve the nonlinear system of equations through Newton's method*)
m = 50;(*number of iterations*)
w = 0.005;(*relaxation parameter*)
init[0] = MapThread[{#1, #2} &, {Join[{ui}, Reverse[Table[((ui - uf)/(b - a)) (i - a) + uf, {i, a + h, b - h, h}]], {uf}], Join[{zi}, Reverse[Table[((zi - zf)/(b - a)) (i - a) + zf, {i, a + h, b - h, h}]], {zf}]}] // Flatten;(*Initial test values*)
For[j = 0, j <= m, j++, residuals = Table[{{resid1}, {resid2}} /. i -> j, {j, 2, n - 1}] // Flatten; 
DFxmat = sparse /. {u[i_] :> init[j][[2 i - 1]], z[i_] :> init[j][[2 i]]}; 
Residvec = Join[{0, 0}, residuals /. {u[i_] :> init[j][[2 i - 1]], z[i_] :> init[j][[2 i]]}, {residbound1, residbound2}/. {u[i_] :> init[j][[2 i - 1]], z[i_] :> init[j][[2 i]]}]; 
init[j + 1] = init[j] + w LinearSolve[DFxmat, -Residvec]] // AbsoluteTiming

(*Residual tolerance of Dirichlet boundary at x=a PLUS interior points*)
ResidTol = Total[((Table[(Abs[resid1]) /. i -> j, {j, 2, n - 1}]) + (Table[(Abs[resid2]) /. i -> j, {j, 2, n - 1}])) /. {u[i_] :> init[j][[2 i - 1]], z[i_] :> init[j][[2 i]]}]/(2 n);
Print["Residual Tolerance = ", ResidTol]
Residual Tolerance = 0.00009350635229

(*Residual tolerance of Dirichlet boundary at x=a PLUS interior points PLUS Neumann boundary at x=b*)
ResidTol = Total[((Table[(Abs[resid1]) /. i -> j, {j, 2, n - 1}]) + (Table[(Abs[resid2]) /. i -> j, {j, 2, n - 1}]) + (Abs[residbound1]) + (Abs[residbound2])) /. {u[i_] :> init[j][[2 i - 1]], z[i_] :> init[j][[2 i]]}]/(2 n);
Print["Residual Tolerance = ", ResidTol]
Residual Tolerance = 0.02602675038

(*Gather data points*)
listu = Take[init[m], {1, -1, 2}];
listz = Take[init[m], {2, -1, 2}];
range = Range[a, b, h];
Listu = MapThread[{#1, #2} &, {range, listu}];
Listz = MapThread[{#1, #2} &, {range, listz}];
uapprox = Interpolation[Listu, InterpolationOrder -> 8, Method -> "Spline"];
zapprox = Interpolation[Listz, InterpolationOrder -> 8, Method -> "Spline"];

(*Plot*)
Plot[{uapprox[x], zapprox[x]}, {x, a, b}, Frame -> True, FrameLabel -> {"x"}, LabelStyle -> Directive[Black, 20], PlotStyle -> {{Green, Thickness[0.005]}, {Blue, Thickness[0.005]}}, PlotRange -> Full, ImageSize -> Large]