How can I solve my system of differential equations?

Question

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["Ersek`RootSearch`"];

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

Answer 1

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

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

---

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.