Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

NDSolve in mathematica drives me crazy and I'm seriously thinking about switching to Matlab because of strange bugs I run into.

NDSolve[{-(2/x) y'[x] - 
     y''[x] - (2 Z/a0 1/x - L (L + 1) - Energy)*y[x] == 0, 
   y[350] == 0, y'[350] == -.0001}, y, {x, .01, 350}] // Flatten
y[1] /. %  

Clearly this should be unsolvable by mathematica because I haven't even defined L or Energy. But for some reason, I'm getting an interpolating function as an output, even after I've cleared all my variables! What on earth is going on!?

I've been having issues trying to find the wavefunctions for the radial part of the schrodinger equation in mathematica, so additionally, if there are any physicists out there, your help showing me a working numeric integration of the radial part of the schrodinger equation would be greatly appreciated.

My full code is:

SolverThingy[n_, L_] :=
   Rydberg = .5;
   En = -1*Rydberg/n^2;
   Z = 1;
   r0 = 2*n*(n + 15);
   rfinal = .5;

   solr = 
    NDSolve[{y''[x] + (L (L + 1)/(2 x^2) + Z/x + En) y[x] == 0, 
       y[350] == 0, y'[350] == -.0001}, y, {x, .01, 350}] // Flatten;
   sol = y /. solr;

   int = Integrate[Evaluate[sol[y]*sol[y]], {y, .01, 350}];
   normd = 1/int*sol[x]
Plot[Evaluate[SolverThingy[10, 1]], {x, .1, 350}, PlotRange -> All}]

For some reason this is NOT the same solution as the analytic solution for the radial equation. And I can't figure out why it doesn't work.

share|improve this question
Running the first code I get the error "NDSolve::ndnum: Encountered non-numerical value for a derivative at x == 350." as expected. Does Context[L] say "Global"? – ssch Oct 17 '13 at 18:44
Plus, as a rule of thumb you should avoid using upper case letters as function names, variables, etc. – Sektor Oct 17 '13 at 18:49
Well..when I copy+paste what I wrote above for the first part, I get an interpolating function and 0.0697003 as a result... I've tried quitting the kernel, exiting mathematica, and resetting variables. – Steven Sagona Oct 17 '13 at 18:52

[Too long for a comment..]

I don't know what result you expect since nothing is shown regarding what would be a reasonable sanity check. Below is fairly minimal code to cover some of the computations of interest.

SolverThingy[n_, L_] := Module[{y, Rydberg = .5, En, Z = 1},
  En = -1*Rydberg/n^2;
  y /. NDSolve[{y''[x] + (L (L + 1)/(2 x^2) + Z/x + En) y[x] == 0, 
      y[350] == 0, y'[350] == -.0001}, y, {x, .01, 350}][[1]]

sol = SolverThingy[10, 1];

int = NIntegrate[sol[y]^2, {y, .01, 350}]

Out[105]= 12.9376

Plot[sol[x], {x, .1, 350}, PlotRange -> All]

enter image description here

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.