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

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?

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

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}

• thank you! I am actually shocked by the grammer style you are using, e.g., when define the function f; how do you know mathematica allows such definitions (with :=  :=)? I want to learn its details very much. If both the left and right boundary conditions have a parameter $x_{\text{left}},x_\text{right}$ to be determined, how do you rewrite it? Why you use MapThread in this case to define the sol1 and sol2? what is the advantage? Why in the FindRoot sentence, there are two brackets with $x_0$s ? The codes work, but I don't know why. I tried to program it similary but failed. – LCFactorization Feb 12 '16 at 12:19
• @LCFactorization f is just a function that creates function. I create this function for conciseness. These 2 lines can be replaced by sol1[x0_?NumericQ] := NDSolveValue[{eqn, bcl, bcr[[1]]}, y, {x, -14, 0}]; sol2[x0_?NumericQ] := NDSolveValue[{eqn, bcl, bcr[[2]]}, y, {x, -14, 0}];, now it's a little easier to understand, right? – xzczd Feb 12 '16 at 13:25
• @LCFactorization See my edit. If you still have difficulty in understanding, evaluate solL[14] alone and check the output. Notice that precisely speaking, $x_0$ is not undetermined, but to be determined. – xzczd Feb 12 '16 at 15:21
• @LCFactorization It's just derivative symbol, when xL gets a numeric value, solL[xl] will evaluate to a InterpolatingFunction. As to the reading material part, have you ever read Leonid Shifrin's excellent book? Here's the unfinished Chinese edition: tieba.baidu.com/p/3230448463 – xzczd Feb 13 '16 at 2:19
• @LCFactorization ?NumericQ is necessary to suppress the warning, if you don't mind the warning, simply use solL[xl_]. There're many examples for the usage of ?NumericQ in this site, you can have a search. – xzczd Feb 13 '16 at 2:44