# NDSolve finding solution when it shouldn't

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

ClearAll["Global*"]
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}];
(*Normalization*)
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.

-
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]
`

-