Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I have a question pertaining to my code below: I require to find the values of the roots E2 that will set the function f1[E2] - f2[E2] to zero. I know how to obtain the roots one at a time: for instance a root of E2 lying between, say 0 and 2. But I would like to know whether it is possible to get Mathematica to generate all the roots that lie within a certain interval at once, say those lying between 0 and 10. I have tried using NSolve and Reduce, but I got nothing, so any help will be most appreciated!

xmin = -6.5; 
xmax = 6.5; 
xmatch = 1.5; 
eq1[x_, x0_, E2_] := {Derivative[2][y][x] + (E2 - x^2) y[x] == 0, 
                      y[x0] == 0, Derivative[1][y][x0] == 1/10^6}; 

y1[x_, E2_] := y[x] /. NDSolve[eq1[x, xmin, E2], y, {x, xmin, xmatch}][[1]]; 
f1[E2_] := Abs[D[y1[x, E2], x]/y1[x, E2] /. x -> xmatch]; 
y2[x_, E2_] := y[x] /. NDSolve[eq1[x, xmax, E2], y, {x, xmax, xmatch}][[1]]; 
f2[E2_] := Abs[D[y2[x, E2], x]/y2[x, E2] /. x -> xmatch]
f[(E2_)?NumericQ] := f1[E2] - f2[E2]; 
(* The following trial doesn't success *)
NSolve[f[(E2_)?NumericQ] == 0 && 0 < E2 < 10, E2]
rts = Reduce[f == 0 && 2.5 < E2 < 10., E2]
c2 = N[E2 /. {ToRules[rts]}]
Print["The value of E2   was=", c2]
share|improve this question
eq1[x_, x0_, E2_] := {(y^\[Prime]\[Prime])[x] + (E2 - x^2) y[x] == 0, y[x0] == 0, Derivative[1][y][x0] == 1/10^6}; is highly suspect. Did you really mean y''[x] rather than (y^\[Prime]\[Prime])[x] –  m_goldberg Dec 10 '13 at 14:53
hi,yes..it should have been y''[x]. –  raj Dec 10 '13 at 15:18
add comment

1 Answer

I suggest you to press Ctrl+Shift+N and Ctrl+Shift+I in your notebook first to change your code into InputForm before you post your code here so it will be cleaner and probably draw more attention. (This time I've done it for you :) )

Then I'd like to point out why your effort leads to nothing. NSolve and Reduce are functions that only work on analytic expression, while InterpolatingFunction generated by NDSolve isn't the case, by the way, the syntax of your equations in NSolve and Reduce are also incorrect, a correct one should be:

f[E2] == 0 && 0 < E2 < 10

Finally, let me fix your code. In fact some posts for similar issue already exists in this site. (Have a look at this question.), but for your question, considering the nature of your equation, we can turn to a even simpler approach, that's all because, your equation can be solved symbolicly:

xmin = -(65/10);
xmax = 65/10;
xmatch = 15/10;
eq[x_, x0_, E2_] = {y''[x] + (E2 - x^2) y[x] == 0, y[x0] == 0, y'[x0] == 1/10^6};

exp = # - #2 & @@ ((Abs[y'[xmatch]/y[xmatch]] /. 
       DSolve[eq[x, #, E2], y, x][[1]]) & /@ {xmin, xmax})
(* The result is a little long so I'd like to omit it here *)

Now what we should do is just solve exp == 0. Unfortunately Solve and NSolve and Reduce are still unavailable, but don't worry, there already exists a robust package called RootSearch for this! After installing it:

rst = RootSearch[exp == 0, {E2, 0, 10}]
{{E2 -> 1.}, {E2 -> 1.76087}, {E2 -> 3.}, 
 {E2 -> 5.}, {E2 -> 5.30345}, {E2 -> 7.}, {E2 -> 9.}}

Some errors generate but the result is reliable:

pts = {Red, Point[{E2, 0} /. rst]}
Plot[exp, {E2, 0, 10}, Epilog -> pts]

enter image description here


For the new problem you mentioned in the comment below, after some trial, I think what you mean is not 2(E2 - x^2), but (2 E2 - x^2) by saying "inserted a factor of 2 in front of (E2-x^2)", right? Then, to resolve it, some modification for my exp is needed. It should be modified like this:

xmin = -(85/10);
xmax = 85/10;
xmatch = 15/10;
eq[x_, x0_, E2_] = {y''[x] + (2 E2 - x^2) y[x] == 0, y[x0] == 0, y'[x0] == 1/10^6};

exp = (# - #2) & @@ ((y'[xmatch]/y[xmatch] /. 
        DSolve[eq[x, #, E2], y, x][[1]]) & /@ {xmin, xmax}) // FunctionExpand
(* The result is a little long so I'd like to omit it here *)

Here, two modifications are made.

First, I removed the Abs function in the exp, though I still don't understand the physics background of this expression very well, my intuition told me that the Abs is just used to avoid very small complexes, but unfortunately it hasn't been used properly so it causes redundant roots, am I right?

Second, which is the most important, is that I've add a FunctionExpand in the definition of exp, this function only changes the form of exp, but with this modification, RootSearch won't lose roots anymore. I thought of adding this after noticed the warning of RootSearch seems to suggest some errors in numerical computation, and a simplification for expressions can sometimes avoid this. Anyway, this time we got the desired result, no fake, no lose:

rst = RootSearch[exp == 0, {E2, 0, 10}]
 {{E2 -> 0.5}, {E2 -> 1.5}, {E2 -> 2.5}, {E2 -> 3.5}, {E2 -> 4.5}, 
  {E2 -> 5.5}, {E2 -> 6.5}, {E2 -> 7.5}, {E2 -> 8.5}, {E2 -> 9.5}}
share|improve this answer
brilliant-tqvm xzczd! never imagined that the equation can be solved symbolically –  raj Jan 20 at 11:45
@raj Glad it helped, and will be more glad if you accept this answer by clicking the tick at the top left corner of my answer :D (If you still feel confused, feel free to continue to ask in the comment :) ) –  xzczd Jan 20 at 12:10
hi.. have a couple of qns . i've noticed that you have used the # & syntax in line number 5 of the code. i don't quite follow . can you kindly explain this to me ? i reran the code for xmin = -8.5, xmax =8.50 , xmatch=1.50 and inserted a factor of 2 in front of (E2-x^2). This is the famous harmonic oscillator problem in quantum physics and the roots are well known =(n+1/2), n=0,1,2,3 etc An error message was generated. i got some of the roots 0.5,1.5, 2.5 ,7.5, 8.5..but root search missed out on 3.5,4.5,5.5. is there another algorithm to get all the roots? cheers raj –  raj Jan 20 at 13:23
@raj For the # & part, it's pure function. You can find very detailed explanation in the document by searching Function in it, also, have a look at this post. For the second question, see my edit. –  xzczd Jan 20 at 14:52
once again, a big thank you, xzczd. you've been a great help. another qn if i may. suppose i have another function, not a nice x^2 function( which gives a pleasant analytical solution in terms of hermite polynomilas ), but something "uglier" like -20/(1 + exp[0.2 x]) without a chance of an analytical solution. How would we be able to generate the roots in this case? My guess would be to use NDsolve and rid the function expand command. I tried that but got no where.any ideas? cheers raj –  raj Jan 21 at 8:00
show 1 more comment

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.