# NDSolve not working for certain parameters

I'm trying to solve the following ODE

myzero = 0.00001;
q = 0.1;
mu = 1;
f[z_] = 1-z^2;
sol = NDSolve[{ f[z] (z y''[z] + y'[z]) == q^2 z y[z],
y[1] == 0, y'[1] == mu}, y, {z, myzero, 1}]


Mathematica is able to solve it for q=0 but for any finite real value it returns the errors Power::infy, Infinity::indet and consequently NDSolve::ndnum (non-numerical value for the derivative). Curious what's the issue. Any help would be appreciated.

• Maple says "Error, (in DEtools/convertsys) unable to convert to an explicit first-order system ". – user64494 Jul 24 '17 at 16:15
• Ignoring error messages doesn't seem wise. (I included a brief synopsis of them.) – Michael E2 Jul 24 '17 at 17:37

The derivatives on the left hand side of your differential equation are multiplied by f[z], which is zero at z=1. This leads to the error messages

Power::infy: Infinite expression 1/0. encountered.
Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered.
NDSolve::ndnum: Encountered non-numerical value for a derivative at z == 1..


because NDSolve seems to divide by f[z] in setting up to solve the equation numerically.

One possible solution is to stop just short of z=1, as you do with z=0 already.

sol = NDSolve[{f[z] (z y''[z] + y'[z]) == q^2 z y[z], y[1 - myzero] == 0,
y'[1 - myzero] == mu}, y, {z, myzero, 1 - myzero}]


seems to work.

• It should be noticed you present a substitute, not a solution for the posed problem (which likely has no solution) @Chris K (which disappears when submitting). – user64494 Jul 24 '17 at 19:26
• Thank you Chris! @user64494 I guess you're right, at z=1 there's singularity, nevertheless for me Chris's solution works :) – Skylar15 Jul 24 '17 at 19:28
• @Skylar15: The question arises: how to apply it to the ill posed problem under consideration? BTW, the ParametricNDSolve command (see reference.wolfram.com/language/ref/ParametricNDSolve.html ) may be applied here – user64494 Jul 24 '17 at 19:40
• @user64494 I suppose I meant "solution" in the colloquial sense of "a means of solving a problem or dealing with a difficult situation", not the mathematical sense ;-) – Chris K Jul 24 '17 at 19:50
• @Chris K: Isn't the use of ParametricNDSolve a better way to this end? – user64494 Jul 25 '17 at 5:26

A solution can be obtained symbolically.

q = 1/10; [z_] = 1 - z^2;


(For some reason q must be rational for FunctionExpand to work properly below.)

y[z] /. Flatten@DSolve[{f[z] (z y''[z] + y'[z]) == q^2 z y[z]}, y[z], z]
(* C[1] Hypergeometric2F1[-(I/20), I/20, 1, z^2] +
C[2] MeijerG[{{}, {1 - I/20, 1 + I/20}}, {{0, 0}, {}}, z^2] *)
sol = FunctionExpand[%, 0 <= z <= 1]
(* C[1] Hypergeometric2F1[-(I/20), I/20, 1, z^2] +
(1 - z^2) C[2] Hypergeometric2F1[1 - I/20, 1 + I/20, 2, 1 - z^2] *)


Now, apply the boundary conditions to sol. (Applying the boundary conditions directly to DSolve fails, probably because MeijerG is an unfriendly function.

sol/. z -> 1
(* {C[1]/(Gamma[1 - I/20] Gamma[1 + I/20]) *)


So, C[1] == 0.

D[sol /. C[1] -> 0, z] /. z -> 1
(* -2*C[2] *)


and C[2] == -1/2.

Plot[sol /. {C[1] -> 0, C[2] -> -1/2}, {z, 0, 1}, AxesLabel -> {z, y},
ImageSize -> Large, LabelStyle -> Directive[12, Bold, Black]]


• sol is defined for $z>0,z<1$, but you substitute $z=1$ in sol ( @bbgodfrey does not appear when submitting). – user64494 Jul 25 '17 at 15:43
• Maple completely confirms your result (The Maple code on demand through Dropbox.). @bbgodfrey does not appear when submitting. +1. – user64494 Jul 25 '17 at 16:02
• @user64494 I presume you are referring to the limits I used with FunctionExpand. I could just as well have used 0 <= z <= 1`, and perhaps should have. (I shall make the change.) Thanks for verifying the solution with Maple and for up-voting the answer. By the way, my user name does not appear when someone writes a comment here, because, as author of the answer, I receive all comments automatically. Best wishes. – bbgodfrey Jul 25 '17 at 16:07
• Here it is as a PDF file dropbox.com/s/47fpdfyz3nfffqk/… – user64494 Jul 25 '17 at 17:26