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]
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}]}
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
