I have a coupled boundary ODE with dependent variables $u=u(x)$ and $z=z(x)$,

$$u'' - \frac{1}{z} \left( -3 + u'^2 (3 - c\; e^{-g u} z^4) - 6 u' z' \right) = 0\tag{1}$$

$$z'' + c\; e^{-g u} z^3 (-3 + 4 u'^2 - 2 u' z') - \frac{3}{z} \left( -1 + u'^2 - 2 u' z' \right) - c^2 e^{-2 g u} z^7 u'^2 - \frac{1}{2} c\; g\; e^{-g u} z^4 u'^2 = 0 \tag{2}$$

where the boundary conditions (B.C.) are fixed on one point $u_f(b)$, $z_f(b)$, while the other is free to move and a derivative condition is imposed,

$$u'_i(a) = \pm \sqrt{\frac{4}{1 - 4 \left(1 - c\; e^{-g u_i(a)} z_i(a)^4 \right)}}\tag{3}$$

$$z'_i(a) = \pm (1 - c\; e^{-g u_i(a)} z_i(a)^4) \sqrt{\frac{4}{1 - 4 \left(1 - c\; e^{-g u_i(a)} z_i(a)^4 \right)}}\tag{4}$$

The domain is given by $x \in [a,b]$. So, the free parameters here are $u_i(a)$ and $z_i(a)$ which we can input "guess" values to start the calculation.

I'm using finite-difference method while employing a relaxation technique to solve the boundary ODE which are eq01 and eq02 in the code.

Parameters: 0 < c < (not too large, maybe at most 10^5) 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 a residual 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 = 1000;
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 Rank 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
```