How to handle a special Neumann-like boundary condition in NDSolve?

Question

How to solve the following nonlinear ODE with two algebraic equations and one boundary condition?

$$y''(x)=\dfrac{2\left((x+15)y'(x)-y(x)\right)\left(y'(x)^2+1\right)}{\left(y(x)^2+x(x+30)+236\right)^2}$$

The Dirichlet boundary condition: $$y(-14)=0$$

The Neumann-like boundary condition with a parameter $x_0$ undetermined:

$$\left\{ \begin{array}{ll} y(x_0)=\sqrt{1-x_0^2} &\\[15pt] y'(x_0)=\dfrac{-x_0}{\sqrt{1-x_0^2}}& \text{where: }-1\lt x_0\lt 0 \\ \end{array} \right.$$

The problem has a physical background in elementary geometric optics therefore the existence and uniqueness of the solution can be easily verified.

Closed form solution might be difficult to obtain. How can NDSolve handle such a boundary-like constraints? Shooting method? Solution domain should therefore be $[-14,x_0]$

In the tutorial document on NDSolve, there is one example on Boundary Value Problems with Parameters, which converts the parameter-contained differential equation into an ODE system by introducing the parameter as a new dependent variable while keeping its first order derivative as zero:

Block[{R = 1}, 
 sol = NDSolve[{f'''[t] - R ((f'[t])^2 - f[t] f''[t]) + R a[t] == 0, 
    a'[t] == 0, f[0] == f'[0] == f'[1] == 0, f[1] == 1}, {f, a}, t];
 Column[{Plot[f[t] /. First[sol], {t, 0, 1}],
   a[0] /. First[sol]}]]

But the method does not solve this problem with a parameter in the boundary conditions instead. How to solve it similarly?

It seems ParametricNDSolve does not handle such problems.

Update

I update this only want to learn the mathematica grammar from @xzczd

If I change the same problem into the following form, i.e. there are parameters in both the left and right boundary conditions of the same nonliear ODE:

left boundary:

$$y(x_l)=-\sqrt{1-(x_l+15)^2},y'(x_l)=\dfrac{x_l+15}{\sqrt{-(x_l+14)(x_l+16)}},\quad -15<x_l<-14$$

right boundary:

$$y(x_r)=\sqrt{9-x_r^2},\quad,y'(x_r)=\dfrac{-x_r}{\sqrt{9-x_r^2}},\quad -3<x_r<0$$

Then how will you rewrite the code per your programming style?

Answer 1

Somewhat surprising that ParametricNDSolve can't handle this type of parameter. Alright, let's go back to the old-fashioned _?NumericQ:

eqn = y''[x] == (2 ((x + 15) y'[x] - y[x]) (y'[x]^2 + 1))/(y[x]^2 + x (x + 30) + 236)^2;

bcl = y[-14] == 0;

bcr = {y[x0] == Sqrt[1 - x0^2], y'[x0] == -x0/Sqrt[1 - x0^2]};

f[sol_, bcr_] := (sol[Pattern[#, _]?NumericQ] := 
                    NDSolveValue[{eqn, bcl, bcr}, y, {x, -14, 0}]) &@x0

MapThread[f, {{sol1, sol2}, bcr}];

FindRoot[sol1[x0][x0] == sol2[x0][x0], {x0, -1/2, -1, 0}]
Plot[sol1[x0][x] /. % // Evaluate, {x, -14, 0}]

{x0 -> -0.0716993}

---

Response to the new added question

The idea is the same:

eqn = y''[x] == (2 ((x + 15) y'[x] - y[x]) (y'[x]^2 + 1))/(y[x]^2 + x (x + 30) + 236)^2;
bcl = {y[xl] == -Sqrt[1 - (xl + 15)^2], y'[xl] == (xl + 15)/Sqrt[-(xl + 14) (xl + 16)]};
bcr = {Sqrt[9 - xr^2], -xr/Sqrt[9 - xr^2]};

(solL[Pattern[#, _]?NumericQ] := NDSolveValue[{eqn, #2}, y, {x, -15, 0}]) &[xl, bcl]

FindRoot[{solL[xl][xr], solL[xl]'[xr]} == bcr, {{xl, -29/2, -15, -14}, {xr, -3/2, -3, 0}}]

{xl -> -14.7652, xr -> -0.816397}