# How can I solve my system of differential equations?

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]


this is the code that i have written to numerically try to solve the DEs to obtain the roots E2. i am able to obtain the roots by plotting "exp" versus E2. however, when i try to use the erzek root finder, an error message is generated

xmin = (-8.5);

xmax = 8.5;

xmatch = 0.50;

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] /.
NDSolve[eq[x, #, E2], y, {x, #, xmatch}][[1]]) & /@ {xmin,
xmax});

Plot[exp, {E2, 0, 10}]

Needs["ErsekRootSearch"];

rst = RootSearch[exp == 0, {E2, 0, 5}](* gives an error message  *)

• 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] Dec 10, 2013 at 14:53
• hi,yes..it should have been y''[x].
– raj
Dec 10, 2013 at 15:18

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]


Update

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 *)


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}}


Update 2

Well, before fixing your new code, I'd like to remind you that, it's not good to completely rewrite your question. (Consider how confused the late-comers be.) Also, after your question being edited for several times, you should have noticed that, you can simply add 4 spaces in front of your codes, or select your codes and press Ctrl+K to format it.

The reason for the error of your new code is simple: Just like Solve, NSolve and Reduce, RootSearch is also a function for analytic equations i.e. it can't handle the InterpolatingFunction generated by NDSolve. As I already mentioned in the comments below, in order to find the roots of an equation containing InterpolatingFunction, the best choice is to we can use ParametricNDSolve or ParametricNDSolveValue together with Ｊ．Ｍ．'s FindAllCrossings:

xmin = (-8.5); xmax = 8.5; xmatch = 0.50;
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] /. y ->
ParametricNDSolveValue[eq[x, #, E2], y, {x, #, xmatch}, E2][E2]) & /@ {xmin, xmax});

ParallelTable[FindAllCrossings[exp, {E2, i - 1, i}, PlotRange -> 1], {i, 10}]

{{0.5}, {1.5}, {2.5}, {3.5}, {4.5}, {5.5}, {6.5}, {7.5}, {8.5}, {9.5}}


Notice that the numeric method isn't that robust (at least when compared to RootSearch) so careful adjustments are made here.

Update 3

OK… after a second look, I noticed that actually RootSearch can handle non-analytic equations under certain conditions. I haven't look into the source code of RootSearch so I can't confirm the exact condition (maybe "the expression should be derivable"?), but it did manage to solve the exp formed by ParametricFunction:

xmin = (-8.5); xmax = 8.5; xmatch = 0.50;
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] /.
y -> ParametricNDSolveValue[eq[x, #, E2], y, {x, #, xmatch},
E2, WorkingPrecision -> 32][E2]) & /@ {xmin, xmax});

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}}


Notice I've added WorkingPrecision -> 32 in exp.

• brilliant-tqvm xzczd! never imagined that the equation can be solved symbolically
– raj
Jan 20, 2014 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 :) ) Jan 20, 2014 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, 2014 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. Jan 20, 2014 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, 2014 at 8:00