If you have a nonlinear second-order boundary value problem where the domain of the problem is $x \in [a,b]$, the boundary conditions imposed are the Robins condition (mixture of the dependent variable z(a) and the derivative z'(b)) at $x=a$ and the Dirichlet condition at $x=b$, we can use the finite difference method to solve this. However, the resources that I've seen only talks about very simple Robins boundary conditions. For example,
\begin{equation} \frac{d^2z}{dx^2} = f(z,z',x) \end{equation}
\begin{equation} c_1 z'(a) + c_2 z(a) = c_3 \end{equation}
We can write these down in terms of finite differences (central for the interior and forward for the left boundary) where the grid points are $i=1 , \dots, n$ and the step size is $h$.
The governing equation is,
\begin{equation} \frac{z_{i + 1} - 2 z_i + z_{i - 1}}{h^2} = f\left(z_i,\left(\frac{z_{i + 1} - z_{i - 1}}{2 h} \right),x_i\right) \quad i=2,\dots,n-1 \end{equation}
Now, the solution $z(x)$ should satisfy the governing equation but it should also satisfy the boundary conditions. So, we need to write down a difference equation at the boundary $x=a$ or at the node $i=1$ that satisfies both the governing equation and the boundary conditions.
Doing a Taylor expansion for $z_2$ and $z_3$ at the point $z_1$,
\begin{align} & z_2 = z_1 + \frac{dz}{dx}\Bigg|_{i=1} h + \frac{d^2z}{dx^2}\Bigg|_{i=1} \frac{h^2}{2!} + \ldots \\ & z_3 = z_1 + \frac{dz}{dx}\Bigg|_{i=1} (2h) + \frac{d^2z}{dx^2}\Bigg|_{i=1} \frac{4h^2}{2!} + \ldots \end{align}
then we have,
\begin{equation} \frac{d^2z}{dx^2}\Bigg|_{i=1} = \frac{-z_3 + 8 z_2 - 7z_1 - 6 h \frac{dz}{dx}\Big|_{i=1}}{2 h^2} \end{equation}
We know from the Robins condition above that,
\begin{equation} \frac{dz}{dx}\Big|_{i=1} = \frac{c_3 - c_2 z_1}{c_1} \end{equation}
and also the governing equation at $i=1$ is,
\begin{align} & \frac{d^2z}{dx^2}\Bigg|_{i=1} = f\left(z_1,\frac{dz}{dx}\Big|_{i=1},x_1\right) \end{align}
Thus, we have
\begin{equation} \frac{-z_3 + 8 z_2 - 7z_1 - 6 h \frac{c_3 - c_2 z_1}{c_1}}{2 h^2} = f\left(z_1,\frac{c_3 - c_2 z_1}{c_1},x_1\right) \end{equation}
Assuming $f$ is not very complicated, i.e. it could be nonlinear but not anything complicated, say a polynomial, then you could collect like terms and maybe have a few nonlinear polynomial terms.
In my situation, the nonlinear second order boundary value differential equation is
\begin{equation} z''(x)-\frac{\frac{1}{100} z(x)^4 \left(2 z'(x)^2+12\right)-600 \left(z'(x)^2+1\right)-\frac{3 z(x)^8}{500000}}{20 z(x) \left(10-\frac{z(x)^4}{1000}\right)}=0 \end{equation}
I have a nonlinear Robins boundary condition at $x=a$. \begin{equation} z'(a) + (d-1) \left(1-\left(\frac{z(a)}{z_h}\right)^{d+1}\right) \sqrt{1+\frac{z'(a)^2}{1-\left(\frac{z(a)}{z_h}\right)^{d+1}}}=0 \end{equation}
where $d=3$, $z_h=10$, and $x\in [10^{-8},10^{-1}]$. Also, the Dirichlet condition at $b=10^{-1}$ is $z(b) = 10^{-3}$.
Discretizing the ODE, we have the governing equation for $i = 2, \ldots, n-1$
\begin{equation} \frac{z_{i + 1} - 2 z_i + z_{i - 1}}{h^2}-\frac{\frac{1}{100} z_i^4 \left(2 \left(\frac{z_{i + 1} - z_{i - 1}}{2 h} \right)^2+12\right)-600 \left(\left(\frac{z_{i + 1} - z_{i - 1}}{2 h} \right)^2+1\right)-\frac{3 z_i^8}{500000}}{20 z_i \left(10-\frac{z_i^4}{1000}\right)}=0 \end{equation}
I have tried discretizing the Robins boundary conditions in the same way as above, i.e. by isolating $z'(a)$, but I'm getting solutions that are doubtful to be correct. So, should I do anything here before doing the finite difference? In particular, I'm thinking the square root is causing some issues. Any guidance on this or suggestion if it does not differ from usual methods? My code is written in Mathematica but I added some comments so that those who are not familiar with it can still kind of follow.
(*Setup the equation*)
Needs["VariationalMethods`"]
f = 1 - (z[x]/zh)^(d + 1);
L = (Sqrt[1 + (z'[x]^2/f)]/z[x]^d) + (d - 1) (z'[x]/z[x]^d);(*Lagrangian*)
eulageq = EulerEquations[L, z[x], x];(*Euler-Lagrange equation*)
s = Solve[eulageq, z''[x]][[1]] // Simplify;
eq = z''[x] - s[[1, 2]] /. {d -> 3, zh -> 10};(*governing equation displayed in post, s[[1,2]] is like the f(z,z',x) in the post*)
bc = z'[x] + (d - 1) (1 - (z[x]/zh)^(d + 1)) Sqrt[1 + (z'[x]^2/(1 - (z[x]/zh)^(d + 1)))];(*Robins condition displayed in post*)
(*Setting up the finite difference and residuals*)
n = 1000;(*Grid points*)
h = (b - a)/(n - 1);(*Step size*)
a = 10^-8;(*a & b are the domain*)
b = 10^-1;
zf = 10^-3;(*zf is the Dirichlet condition at x=b*)
zp = -Sqrt[((d - 1)^2 (1 - (z[x]/zh)^(d + 1))^2)/(1 - (d - 1)^2 (1 - (z[x]/zh)^(d + 1)))] /. {d -> 3, zh -> 10, z[x] -> z[1]};(*zp is the z'[a] solved from the boundary conditions bc*)
rule = 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}];(*finite difference rule*)
eqs = Table[{eq} /. rule[[i]], {i, Length[rule]}];(*substitute the finite difference to the governing equation*)
residual = h^2 eqs // Flatten;(*residual of the governing equation*)
residbound = (-z[3] + 8 z[2] - 7 z[1] - 6 zp h) - 2 h^2 s[[1, 2]] /. {d -> 3, zh -> 10, z[x] -> z[1], z'[x] -> zp};(*residual of the Robins condition, s[[1,2]] is like the f(z,z',x) in the post*)
(*Setup the sparse matrix*)
For[i = 2, i <= n - 1, i++, jac[i] = D[residual[[i - 1]], {{z[i - 1], z[i], z[i + 1]}}]]
DFx = Table[jac[i], {i, 2, n - 1}];
ShiftMatrix[mat_, shift_] := Reverse@PadLeft@MapThread[PadLeft[#1, Length[mat] + #2, 0, #2] &, {Reverse[mat], shift}]
jacbound = D[residbound, {{z[1], z[2], z[3]}}];
sparseresidual = ShiftMatrix[DFx, Table[i, {i, 0, n - 3}]][[All, n - 4 ;;]];
sparse = Join[{Join[jacbound, ConstantArray[0, n - 3]]}, sparseresidual, {Join[ConstantArray[0, n - 1], {1}]}];
m = 90;(*Number of iteration*)
zi = 9.5;(*Initial test point for z[1]*)
z0[0] = Join[{zi}, Reverse[Table[((zi - zf)/(b - a)) (i - a) + zf, {i, a + h, b - h, h}]], {zf}];(*Initial test points for the z[i]'s*)
(*Newton's method*)
For[j = 0, j <= m, j++, residuals = h^2 eqs;
DFxmat = sparse /. z[i_] :> z0[j][[I]];
Residvec = {residbound /. z[i_] :> z0[j][[i]], residuals /. z[i_] :> z0[j][[i]], 0} // Flatten;
z0[j + 1] = z0[j] + 0.22 LinearSolve[N[DFxmat], N[-Residvec]] // Flatten] // AbsoluteTiming
(*Residual error*)
ResidTol = Total[Table[Abs[residuals[[i]]], {i, 1, n - 2}] /. z[i_] :> z0[j][[i]] // Flatten]/n;
Print["Residual Tolerance = ", ResidTol]
Residual Tolerance = 7.090653622*10^-14
range = Range[a, b, h];
list = MapThread[{#1, #2} &, {range, z0[j]}];
zapprox = Interpolation[list, InterpolationOrder -> 10, Method -> "Spline"];
zapprox[a]
9.999999722
zapprox'[a]
0.01664915436
zapprox''[a]
-498.9821196
(*Check if Robins boundary is satisfied, must be equal to 0*)
bc /. {d -> 3, zh -> 10, z[x] -> zapprox[a], z'[x] -> zapprox'[a]}
0.01666025571
(*Check if governing equation is satisfied, must be equal to 0*)
zapprox''[a] - s[[1, 2]] /. {d -> 3, zh -> 10, z[x] -> zapprox[a], z'[x] -> zapprox'[a]}
-0.003923034814
(*Using zapprox[a] in the Robins condition bc as a consistency check for zapprox'[a]*)
sol = (z' + (d - 1) (1 - (z[x]/zh)^(d + 1)) Sqrt[1 + (z'^2/(1 - (z[x]/zh)^(d + 1)))]) /. {d -> 3, zh -> 10, z[x] -> zapprox[a]};
Solve[sol == 0, z']
z' -> -2.222100815*10^-7
The governing equation and boundary equation at $x=a$ is not 0 (as you see all the terms are on one side so that when I calculate the result should be 0 or at least machine precision). Also, the resulting zapprox'(a) is not consistent with the result given by the Robins condition bc if you plug in zapprox(a).