I have a nonlinear coupled second order differential equation which is going to be generated by solving the Euler Lagrange equation for which those differential equations correspond to eq01 and eq02. I'm going to use finite-difference method to solve these equations, i.e., I'll set up the sparse matrix which when solved will provide the correction needed for the initial guess (indicated by init[0]) of the solution to the differential equations. After obtaining the sparse matrix, I solve it using Newton's method. The code below is divided into these three steps for clarity.
The boundary conditions of the differential equations are ui,uf,zi,zf. The constraint is ui < uf while zf = 10^-1 is fixed. My plan is to run different values of ui and zi for a fixed uf (the difference of uf from zf is I'll eventually pick another value for uf then run different values of ui and zi again but zf is completely fixed for all these). If I choose say, uf = 7 then ui = 5 and zi = 8, everything works fine.
However, when I choose ui or zi to be far away (or both are far away) from uf I encounter an error from LinearSolve. For example, ui = 2 (a bit far from uf = 7) and zi = 8 OR ui = 6 and zi = 15 (a bit far from uf = 7) OR ui = 2 and zi = 15 (both are a bit far from uf = 7). I understand that I can't just take any arbitrary values for these parameters, but at least I want to have some considerable amount of range to run these parameters, e.g. maybe zi=1 to zi=20.
(*Equation setup*) Clear["Global`*"] Needs["VariationalMethods`"] u0 = 20; 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) + (d - 1) (z'[x]/z[x]^d);(*Lagrangian*) 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]]; (*Finite difference method: Sparse matrix setup*) n = 200; (*grid points*) h = (b - a)/(n - 1); (*step size*) a = 10^-8; (*a & b are the domain*) b = 10^-1; ui = 2; (*ui, uf, zi, zf are the boundary conditions*) uf = 7; zi = 15; zf = 10^-1; 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; 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]]}; 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 - 2], {1, 0}]}, {Join[ConstantArray[0, 2 n - 1], {1}]}]; (*Solving the nonlinear system of equations through Newton's method*) m = 30; (*Number of iterations*) 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; 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]]}, {0, 0}]; init[j + 1] = init[j] + 0.7 LinearSolve[N[DFxmat, 30], N[-Residvec, 30]] // Flatten] // AbsoluteTiming finalsolution = init[j] LinearSolve::luc: Result for LinearSolve of badly conditioned matrix {{1.,0.,0.,0.,0.,0.,0.,0.,0.,0.,<<390>>},{0.,1.,0.,0.,0.,0.,0.,0.,0.,0.,<<390>>},<<7>>,{0.,0.,0.,0.,0.,0.,132.1319891,0.9584233208,0.01115998895,-2.012128375,<<390>>},<<390>>} may contain significant numerical errors. LinearSolve::luc: Result for LinearSolve of badly conditioned matrix {{1.,0.,0.,0.,0.,0.,0.,0.,0.,0.,<<390>>},{0.,1.,0.,0.,0.,0.,0.,0.,0.,0.,<<390>>},<<7>>,{0.,0.,0.,0.,0.,0.,105.0729581,0.9683577215,0.007541112841,-2.008287291,<<390>>},<<390>>} may contain significant numerical errors. LinearSolve::luc: Result for LinearSolve of badly conditioned matrix {{1.,0.,0.,0.,0.,0.,0.,0.,0.,0.,<<390>>},{0.,1.,0.,0.,0.,0.,0.,0.,0.,0.,<<390>>},<<7>>,{0.,0.,0.,0.,0.,0.,-1.083697194*10^9,189.6354636,1.494207434*10^6,-249300.3218,<<390>>},<<390>>} may contain significant numerical errors. General::stop: Further output of LinearSolve::luc will be suppressed during this calculation. ``` 
