I have a coupled boundary ODE represented by eq01 and eq02. I'm using finite-difference method while employing a relaxation technique to solve the boundary ODE. The boundary conditions (B.C.) are fixed on one point represented by uf and zf, while the other is free to move and a derivative condition is imposed represented by up and zp. The free moving point is represented by ui and zi, they are "guess" values we need to input to start the calculation.

Parameters: 0 < c < (not too large, maybe at most 10^4) and g is a damping factor which as can be seen is in the exponent; a b zf uf are fixed; n is the number of grid points while h is the grid size; zi ui can be changed so as to lead to convergence.

Procedure (I think this is not necessary to know to understand the issue): I have built a sparse matrix where the free B.C. is on the upper left and the fixed B.C. is on the lower right. resid indicates the residual of the finite-difference not on the boudary while residbound indicates the residual of the boundary. The sparse matrix sparse was constructed, then a relaxation technique was employed with iteration m and relaxation parameter 0 < w < 1; w can be changed for faster convergence and m can be increased for smaller error. In addition, initial values represented by init[0] were injected (it includes the aforementioned ui zi). The matrix DFxmat and the vector Residvec are the sparse matrix and residual vector where inputs were injected.

The issue: Given the matrix DFxmat and the vector Residvec, I used LinearSolve to compute the solution. However, it gives an error saying that DFxmat is badly conditioned. One check that can be done if it is rank deficient, probably due to very small numbers is to calculate MatrixRank = 200 which is the same as the size of the matrix so it is not deficient. However, the condition number gives a really large value ConditionNumber = 10^11. One thing I'm thinking about is maybe it is due to a large difference in the entries of the matrix where some entries could be of order $10^0$ while others are $10^6$. Changing w to some other number may not give a warning error but the ConditionNumber is still big, in addition it still does not give a satisfactory solution (see Addendum). Overall, I'm not sure what is contributing to this issue and I don't know how to resolve this.

Addendum: The solution that I want should give an error of order $10^{-10}$ or at least close. However, the issue is giving an error of $\sim 0.26$. I know that by increasing the grid points n a better approximation can be achieved, but for the purposes of the issue I just set it to n=100.

(*Two Coupled ODE Setup*)
Clear["Global`*"]
Needs["VariationalMethods`"]
c = 10;
g = 3;
m = c Exp[-g u[x]];
f = 1 - m z[x]^(d + 1);
L = Sqrt[-f u'[x]^2 + 2 u'[x] z'[x] + 1]/z[x]^d;
eq1 = EulerEquations[L, u[x], x];
eq2 = EulerEquations[L, z[x], x];
s = Solve[{eq1, eq2} /. d -> 3, {u''[x], z''[x]}][[1]] // Simplify;
eq01 = u''[x] - s[[1, 2]];
eq02 = z''[x] - s[[2, 2]];

(*Parameters*)
n = 100;
h = (b - a)/(n - 1);
a = 0;
b = 10^-2;
zf = 10^-2;
uf = 0.03;
ui = 0.1;
zi = 0.26;
up = Sqrt[Sqrt[((d - 1)^2/(1 - (d - 1)^2 (1 - c Exp[-g u[1]] (z[1])^(d + 1))))^2]] /. d -> 3;
zp = -(1 - c Exp[-g u[1]] (z[1])^(d + 1)) Sqrt[Sqrt[((d - 1)^2/(1 - (d - 1)^2 (1 - c Exp[-g u[1]] (z[1])^(d + 1))))^2]] /. d -> 3;

(*Sparse Matrix Setup*)
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]};
resid1 = h^2 eq01 /. rule;
resid2 = h^2 eq02 /. rule;
residbound1 = ((-u[3] + 8 u[2] - 7 u[1] - 6 up h) - 2 h^2 s[[1, 2]]) /. {d -> 3, u[x] -> u[1], z[x] -> z[1], u'[x] -> up, z'[x] -> zp};
residbound2 = ((-z[3] + 8 z[2] - 7 z[1] - 6 zp h) - 2 h^2 s[[2, 2]]) /. {d -> 3, u[x] -> u[1], z[x] -> z[1], u'[x] -> up, z'[x] -> zp};
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]]};
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[1], z[1], u[2], z[2], u[3], z[3]}}];
parresidbound2 = D[residbound2, {{u[1], z[1], u[2], z[2], u[3], z[3]}}];
mat = {parresid1, parresid2};
sparseresidual = Normal[SparseArray[Table[Band[{2 (i - 2) + 1, 2 (i - 2) + 1}] -> {mat}, {i, 2, n - 1}]]];
sparse = Join[{Join[parresidbound1, ConstantArray[0, 2 n - 6]]}, {Join[parresidbound2, ConstantArray[0, 2 n - 6]]}, sparseresidual, {Join[ConstantArray[0, 2 n - 2], {1, 0}]}, {Join[ConstantArray[0, 2 n - 1], {1}]}];

(*Relaxation Method*)
m = 10;
w = 0.3;
init[0] = MapThread[{#1, #2} &, {Join[{ui}, Reverse[Table[((ui - uf)/(b - a)) (i - a) + uf, {i, a + h, b - h, h}]], {uf}], Table[zi (1 - i^2), {i, 0, Sqrt[1 - zf/zi], Sqrt[1 - zf/zi]/(n - 1)}]}] // Flatten;
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[{residbound1, residbound2}, residuals, {0, 0}] /. {u[i_] :> init[j][[2 i - 1]], z[i_] :> init[j][[2 i]]} // N; init[j + 1] = Chop[N[init[j], 30]] - w LinearSolve[Chop[N[DFxmat, 30]], Chop[N[Residvec, 30]]]//Flatten]

LinearSolve::luc: Result for LinearSolve of badly conditioned matrix {{-6.999999864,0.0001784600311,8.,0.,-1.,0.,0.,0.,0.,0.,<<190>>},{6.802416307*10^-8,-6.999818669,0.,8.,0.,-1.,0.,0.,0.,0.,<<190>>},<<8>>,<<190>>} may contain significant numerical errors.
LinearSolve::luc: Result for LinearSolve of badly conditioned matrix {{-6.999999999,0.002946668739,8.,0.,-1.,0.,0.,0.,0.,0.,<<190>>},{2.876690113*10^-10,-6.997053283,0.,8.,0.,-1.,0.,0.,0.,0.,<<190>>},<<8>>,<<190>>} may contain significant numerical errors.
LinearSolve::luc: Result for LinearSolve of badly conditioned matrix {{-6.999999641,0.0001058640351,8.,0.,-1.,0.,0.,0.,0.,0.,<<190>>},{1.795481327*10^-7,-6.999888067,0.,8.,0.,-1.,0.,0.,0.,0.,<<190>>},<<8>>,<<190>>} may contain significant numerical errors.

(*Matrix Rand and Condition Number*)
MatrixRank[DFxmat]
LinearAlgebra`Private`MatrixConditionNumber[DFxmat, Norm -> 1]
200
1.403654866*10^11

(*Residual Tolerance - Error of the solution*)
ResidTol = Total[Flatten[{Abs[residbound1] + Abs[residbound2] + Table[(Abs[resid1] + 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.2595533347
```