First, there is a mistake in the sign `k`. And secondly, for such tasks we use a special method that expands the possibilities of the shooting method. I will demonstrate the working code A1 = 200; A2 = 186*10^5; A3 = 1/5000; a = Rationalize[A2*A3^3]; H = 1; r0 = 10^-5; f[z_, r_] := A1*(z + H) + a*(1/(Sqrt[z^2 + r^2] - 1))^3; k = -z''[r]/(z'[r]^2 + 1)^(3/2) - z'[r]/(r Sqrt[z'[r]^2 + 1]); sol = ParametricNDSolveValue[{k == f[z[r], r], z'[r0] == 0, z[r0] == z0}, z, {r, r0, 1}, {z0}]; Plot[Evaluate[Table[sol[z0][r], {z0, -1.15, -1.05, .01}]], {r, r0, 1}, PlotRange -> All] [![fig1][1]][1] Find the parameter that satisfies the boundary condition on the right border FindRoot[sol[z0][1] == -H, {z0, -1.05}] Out[]= {z0 -> -1.03176} {Plot[sol[-1.03176][r], {r, r0, 1}],Plot[{sol[-1.03176][r], -Sqrt[1 - r^2]}, {r, r0, 1}]} [![fig2][2]][2] The author insists that the correct model corresponds to his code. We give a solution for this case. The solution method does not change. A1 = 200; A2 = 186*10^5; A3 = 1/5000; a = Rationalize[A2*A3^3]; H = 1; r0 = 10^-5; f[z_, r_] := A1*(z + H) + a*(1/(Sqrt[z^2 + r^2] - 1))^3; k = z''[r]/(z'[r]^2 + 1)^(3/2) + z'[r]/(r Sqrt[z'[r]^2 + 1]); sol = ParametricNDSolveValue[{k == f[z[r], r], z'[1] == z1, z[1] == -H}, z, {r, r0, 1}, {z1}, WorkingPrecision -> 30]; {Plot[Evaluate[ Table[sol[z1][ r], {z1, .00024816068, .00024816069, .000000000001}]], {r, r0, 1}, PlotRange -> All], Plot[Evaluate[ Table[sol[z1][ r], {z1, .00024816068, .00024816069, .000000000001}]], {r, r0, .01}, PlotRange -> All]} // Quiet [![fig3][3]][3] [1]: https://i.sstatic.net/bddoj.png [2]: https://i.sstatic.net/Jn5uP.png [3]: https://i.sstatic.net/UkS7L.png