Jacobian General::ivar issue

Question

I have a differential equation labelled by eq0 and I want to use finite-difference method to solve it. First, I have to set up the system of equations eqn[i] through a For function. Next, I have to calculate the Jacobian of the coefficient matrix of the system of equations for Newton's method. The parameter n labels the number of equations (or the grid points).

When I calculate the Jacobian up to n=8 there is no problem, but beyond that, say n=10 I encountered an General::ivar issue. I tried calculating the Jacobian using the derivative D function but I encountered the same issue so I guess maybe there is something wrong with how I wrote the code. Any hints?

ClearAll["Global`*"]
Needs["VariationalMethods`"]
f = 1 - (z[x]/zh)^(d + 1);
L = Sqrt[1 + (z'[x]^2/f)]/z[x]^d;(*Lagrangian*)
eulageq = EulerEquations[L, z[x], x];(*Euler-Lagrange equation*)
s = Solve[eulageq, z''[x]][[1]] // Simplify;(*2nd order EOM*)
eq0 = z''[x] - s[[1, 2]] /. {d -> 3, zh -> 10};
(Setting up the nonlinear system of equations)
a = 0;
b = 1;
n = 10;
h = (b - a)/(n + 1);
alpha = 95/10;
beta = 10^-3;
z[0] = alpha;
z[n + 1] = beta;
For[i = 1, i <= n, i++, eqn[i] = Simplify[Collect[eq0 /. {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]}, z[i]]]; Print["eqn[", i, "] = ", eqn[i]]]
j = 0;
x[0] = Table[(1 - i) alpha, {i, 1/10, 90/100, (90/100 - 10/100)/(n - 1)}]; 
xr[j] = MapThread[#1 -> #2 &, {Array[z, Length[Table[i, {i, 1, n}]]], x[j]}]; 
DFx = ResourceFunction["JacobianMatrix"][Table[eqn[i], {i, 1, n}], Table[z[i], {i, 1, n}]] /. xr[j]//N
General::ivar: 1/1000 is not a valid variable.
General::ivar: 171/20 is not a valid variable.
General::ivar: 8.55` is not a valid variable.
General::stop: Further output of General::ivar will be suppressed during this calculation.
JacobianMatrix [{73.26472605,46.03916874,45.22798136,47.78091888,53.1565758,61.92990427,75.89993511,419.515069,262284.3839,-226.5003625},{8.55,7.705555556,6.861111111,6.016666667,5.172222222,4.327777778,3.483333333,2.638888889,1.794444444,0.95}]
(Jacobian using D function)
DFx1=D[Table[eqn[i],{i,1,n}],{Table[z[i],{i,1,n}]}]/.xr[j]//N
During evaluation of In[222]:= General::ivar: 1/1000 is not a valid variable.
During evaluation of In[222]:= General::ivar: 171/20 is not a valid variable.
During evaluation of In[222]:= General::ivar: 8.55is not a valid variable. During evaluation of In[222]:= General::stop: Further output of General::ivar will be suppressed during this calculation. \!\( \*SubscriptBox[\(\[PartialD]\), \({{8.55, 7.705555555555556, 6.861111111111111, 6.016666666666667, 5.1722222222222225, 4.3277777777777775, 3.4833333333333334, 2.638888888888889, 1.7944444444444445, 0.95}}\)]\({73.26472605314348, 46.039168736484946, 45.22798135742989, 47.78091887685099, 53.156575796563885, 61.92990426583507, 75.89993511371034, 419.51506904343574, 262284.38392592594, (-226.5003624756579`)}))

Answer 1

With a small correction we can compute solution using FindRoot as follows

ClearAll["Global`*"]
Needs["VariationalMethods`"]
f = 1 - (z[x]/zh)^(d + 1);
L = Sqrt[1 + (z'[x]^2/f)]/z[x]^d;(*Lagrangian*)
eulageq = EulerEquations[L, z[x], x];(*Euler-Lagrange equation*)
s = Solve[eulageq, z''[x]][[1]] // Simplify;(*2nd order EOM*)
eq0 = z''[x] - s[[1, 2]] /. {d -> 3, zh -> 10};
(Setting up the nonlinear system of equations)
a = 0;
b = 1;
n = 250;
h = (b - a)/(n + 1);
alpha = 95/10;
beta = 10^-3; Z = Array[z, {n + 2}];
z[1] = alpha;
z[n + 2] = beta; rul = 
 Table[{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]}, {i, 2, 
   n + 1}];
eqs = Table[eq0 /. rul[[i]], {i, Length[rul]}];
DFx = ResourceFunction["JacobianMatrix"][eqs, 
   Table[z[i], {i, 2, n + 1}]];
x[0] = Table[(1 - i) alpha, {i, 1/10, 90/100, (90/100 - 10/100)/(n)}];
sol = FindRoot[eqs, Table[{z[i], x[0][[i]]}, {i, 2, n + 1}], 
  Jacobian -> DFx, MaxIterations -> 1000, WorkingPrecision -> 30]

Visualization

ListLinePlot[Z /. sol, PlotRange -> All]

Using code from my answer here we can compute solution for inverse function $x(z)$ as follows

Clear[f, x, z]
f = 1 - (z/zh)^(d + 1);
L = (Sqrt[1 + (x'[z]^(-2)/f)]/z^d) x'[z];(Lagrangian)
eulageq = EulerEquations[L, x[z], z];(Euler-Lagrange equation)
s = Solve[eulageq, x''[z]][[1]] // Simplify;
eq0 = x''[z] - s[[1, 2]] == 0;
eq1 = eq0 /. {zh -> 10, d -> 3};
zs = 95/10; z0 = 10^-3; sol2 = 
 NDSolveValue[{eq1, x[zs] == 0, x[z0] == 1}, x, z]

Finally we compare sol and sol2 in one plot

Show[Plot[sol2[z], {z, 10^-3, zs}, PlotRange -> All, 
  AxesLabel -> {"z", "x"}], 
 ListPlot[Table[{z[i] /. sol, xx[[i]]}, {i, n + 2}], 
  PlotStyle -> Red]]