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
```