NDsolve complicated function

I'm trying to solve a complicated second-order ODE with NDSolve. I've never used NDSolve before. I've looked a lot of forums and tried some things, but I haven't been able to get the ODE's solution. My code is

(* Constants *)
α = 0.04;
σ = 0.1;
A = 0.1;
θ = 10;
δ = 0;

(* Function *)
ϕ[i_] := i - θ/2 i^2 - δ

(* Calculated Boundary p *)
i = 0.5 (A + 1/θ - Sqrt[(A - 1/θ)^2 + 4 α/θ]);
p = (A - i) Exp[(ϕ[i] - σ^2/2)/α];

(* Symbolic Functions for ODE *)
L0 = (1 - z) (1 - z y'[z]/y[z]);
L1 = z (1 + (1 - z) y'[z]/y[z]);
c = 0.5 (A - 1/θ) + 0.5 Sqrt[(A - 1/θ)^2 + 4 α/θ ((1 - z)^2/L0 + z^2/L1)];
i0 = (1 - α/c (1 - z y'[z]/y[z])^(-1))/θ;
i1 = (1 - α/c (1 + (1 - z) y'[z]/y[z])^(-1))/θ;
m = z^2 (1 - z)^2 y''[z]/y[z];

(* Solve ODE *)
NDSolve[{α Log[c/y[z]] + ϕ[i0] L0 + ϕ[i1] L1 - σ^2/2 (L0^2 + L1^2) + σ^2 m == 0, y[0] == p, y[1] == p}, y, {z, 0, 1}]


All I get is "Infinite Expression" and "Indeterminate Expression". Is this function too difficult for Mathematica?

Thank you, Jane

Edit: The output is just NDSolve[...] with all the input just spit out again.

• Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. Sep 30 '15 at 5:21
• The ODE seems to be singular at both end points. Sep 30 '15 at 5:59
• Side note: Comments in Mathematica are placed between parentheses + asterisks. What you had before, words followed a semicolon, was actually a computation representing the multiplication of the words (and in general would have to be a syntactically correct expression, which would be computed, however long it takes). Sep 30 '15 at 12:55

2 Answers

Some progress can be made as follows in solving this equation. First, replace all \$MachinePrecision coefficients by exact numbers.

α = 4/100; σ = 1/10; A = 1/10; θ = 10; δ = 0;

ϕ[i_] := i - θ/2 i^2 - δ

i = (A + 1/θ - Sqrt[(A - 1/θ)^2 + 4 α/θ])/2;
p = (A - i) Exp[(ϕ[i] - σ^2/2)/α];

L0 = (1 - z) (1 - z y'[z]/y[z]);
L1 = z (1 + (1 - z) y'[z]/y[z]);
c =  (A - 1/θ)/2 + Sqrt[(A - 1/θ)^2 + 4 α/θ ((1 - z)^2/L0 + z^2/L1)]/2;
i0 = (1 - α/c (1 - z y'[z]/y[z])^(-1))/θ;
i1 = (1 - α/c (1 + (1 - z) y'[z]/y[z])^(-1))/θ;
m = z^2 (1 - z)^2 y''[z]/y[z];


Then, as suggested by Alexei Boulbitch, Simplify the ODE

eq = Simplify[y[z]^2
(α Log[c/y[z]] + ϕ[i0] L0 + ϕ[i1] L1 - σ^2/2 (L0^2 + L1^2) + σ^2 m)] == 0
(* ((5 + 2*z - 2*z^2 + 8*Log[Sqrt[(y[z]*(y[z] + (1 - 2*z)*Derivative[1][y][z]))/
((y[z] - (-1 + z)*Derivative[1][y][z])*(y[z] - z*Derivative[1][y][z]))]/
(5*Sqrt[10]*y[z])])*y[z]^2 - 2*(-1 + z)^2*z^2*Derivative[1][y][z]^2 +
2*(-1 + z)*z*y[z]*((-1 + 2*z)*Derivative[1][y][z] +
(-1 + z)*z*Derivative[2][y][z]))/200 == 0 *)


As I noted in my comment last night, the equation is singular at both end points. Nonetheless, it seems likely that the end points can be handled by expansions in z (as, for instance, Bessel's equation can be integrated numerically near the origin). Here, to obtain a first cut at a solution, we simply move in from these endpoints by a small amount.

s = First@NDSolve[{eq, y[1/100] == p, y[99/100] == p}, y, {z, 1/100, 99/100},
Method -> {"Shooting", "StartingInitialConditions" -> {y[1/100] == p,
y'[1/100] == Rationalize[0.1536914, 10^-10]}}, WorkingPrecision -> 30];
Plot[y[z] /. s, {z, 0.01, .99}, AxesLabel -> {y, z}]


• Wow. That's fantastic help. I was aware of the singularity at the end points and I had tried perturbing them by epsilon but I couldn't proceed. Thank you. Sep 30 '15 at 16:05

If you evaluate your equations together with its simplification:

eq = (α Log[c/y[z]] + ϕ[i0] L0 + ϕ[
i1] L1 - σ^2/2 (L0^2 + L1^2) + σ^2 m // Expand //
Simplify) == 0


you will find there such kind of terms:

The denominator under the logarithm and the square root turns into zero at z->1; y[1]->y'[1]->p:

(y[z] - (-1 + z) Derivative[1][y][z]) (y[z] -
z Derivative[1][y][z]) /. {y[z] -> p, y'[z] -> p, z -> 1}

(*  0. *)


So, the problem is in mathematics, rather than in Mathematica.

Have fun!

• Isn't the problem rather the factors z^2 (1 - z)^2 multiplying y''[z], which bbgodfrey hints at in a comment? Sep 30 '15 at 13:09
• @MichaelE2 Indeed, y'[1] -> p is not correct. The slope is negative near the upper endpoint. Sep 30 '15 at 13:40