# NDSolve : Second order ODE with variable coefficients

I have a second order ODE of the form :

u''[x]+ p[x]u'[x]+ q[x] u[x]==0 ;


with,

p[x]= -(((1.15+ 4.3 x - 6 x^2))/(-0.03333333333333333 + 1.15 x + 2.15 x^2 - 2 x^3))
q[x]= ((6 -(0.06666666666666667/x^2) +(1.15/x) + 2 x))/(-0.03333333333333333+ 1.15 x + 2.15 x^2 - 2 x^3)


With Boundary condition

u[0]==0,u'[0]==0


When I solved it using NDSolve, I got the following errors :

Power::infy: Infinite expression 1/0.^2 encountered. >>
Power::infy: Infinite expression 1/0. encountered. >>
Infinity::indet: Indeterminate expression 6.+ComplexInfinity+ComplexInfinity
encountered. >>
NDSolve::ndnum: Encountered non-numerical value for a derivative at x == 0.. >>


I have used the following code:

Clear[u, t, x, eq, sol, v];
L = 6; m = -0.1; n = 1.15; k = 1.15;
p = (-2 x^3 + (x^2) (1 + k) + n x + m/3);
eq = -((L + (2 x + (n/x) + (2 m)/(3 x^2)) )/p) u[x]
-D[p, x]/p D[u[x], x] +   D[D[u[x], x], x];
eqn = {eq == 0, u[0] == 0, D[u[x],x]/.x->0 == 0};
sol = NDSolve[eqn, u, {x, 0, 1}]

• Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour! 3) 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. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! Oct 17 '15 at 15:39
• You can format inline code and code blocks by selecting the code and clicking the {} button above the edit window. The edit window help button ? is also useful for learning how to format your questions and answers. You may also find this this meta Q&A helpful Oct 17 '15 at 15:39
• Isn't the solution $u(t) \equiv 0$? Oct 17 '15 at 15:47
• No, actually the solution may come in the range of the order of 10^-3 or so. Here is the code i have used : Clear[u, t, x, eq, sol, v] L = 6; m = -0.1; n = 1.15; k = 1.15; p = (-2 x^3 + (x^2) (1 + k) + n x + m/3); eq = -((L + (2 x + (n/x) + (2 m)/(3 x^2)) )/p) u[x] - D[p, x]/p D[u[x], x] + D[D[u[x], x], x]; eqn = {eq == 0, u[0] == 0, u'[0] == 0}; sol = NDSolve[eqn, u, {x, 0, 1}] Oct 17 '15 at 15:50
• Please edit your post with the code, properly formatted, rather than put it in a comment. Comments are by design transitory on stack exchange, so it's better to have all the info in the post itself. Oct 17 '15 at 16:02

Using the original code

eq = u''[x] + p[x] u'[x] + q[x] u[x];

p[x] = -(((1.15 + 4.3 x - 6 x^2))/(-0.03333333333333333 + 1.15 x +
2.15 x^2 - 2 x^3)) // Rationalize;
q[x] = ((6 - (0.06666666666666667/x^2) + (1.15/x) +
2 x))/(-0.03333333333333333 + 1.15 x + 2.15 x^2 - 2 x^3) // Rationalize;


we can use frobeniusNDSolve, which applies Frobenius' method, from my answer to Attempting to use NDSolve to plot harmonic oscillator solutions, to find two independent solutions.

usol = frobeniusNDSolve[eq, u, {x, 0, 1}]


Note that the indicial equation has two conjugate complex roots ${1\over2}\big(1\pm i\sqrt{7}\big)$. We can get two independent real solutions by taking the real and imaginary parts of the Frobenius solution corresponding to one of the roots.

Plot[(Through[{Re, Im}@u[x]] /. First@usol) // Evaluate, {x, 0, 1}]
Plot[(Through[{Re, Im}@u[x]] /. First@usol) // Evaluate, {x, 0, 0.01}]


While the value of u[0] is 0 in both cases, the derivative u'[x] approaches infinity (and is roughly asymptotic to 1/Sqrt[x]).

{Through[{Re, Im}@u[0]], Through[{Re, Im}[u'[1*^-12]]]} /. First@usol
(*  {{0., 0.}, {-999941., -1.00006*10^6}}  *)


Consequently, the only solution satisfying the initial condition u[0] == u'[0] == 0` is the trivial solution.

• @Pisces You're welcome. Oct 21 '15 at 12:32
• I've checked the above code with the help of frobeniusNDSolve. But it shows error. Oct 21 '15 at 15:27
• @Pisces I'll take a look when I get a chance. I updated the code earlier, but maybe I forgot to include something. Sorry about the inconvenience. Oct 21 '15 at 15:47
• Yea, Sure and no inconvenience at all. You have already helped me a lot. Thanks so much for spending your time Oct 21 '15 at 15:51
• @Pisces I checked with a fresh kernel and the code worked without trouble. Some possibilities: Maybe you have a lurking definition that interferes with my code? Perhaps you made a change or tried a different example? What version of Mathematica are you using? What errors do you get? Oct 22 '15 at 10:58