Bug introduced in Version 11 or earlier and persisting through 12.1. Reported to Wolfram Technical Support as CASE:4532301.

I am considering the eigenvalue problem associated with the double-well harmonic oscillator. Using DEigensystem,

DEigensystem[-1/2 y''[x] + (-x^2/2 + x^4/4) y[x], y[x], {x, -∞, ∞}, 2]

gives eigenvalues that depend on $x$:

{Sqrt[-2 + x^2]/(2*Sqrt[2]), (3*Sqrt[-2 + x^2])/(2*Sqrt[2])

Maybe Mathematica 11 cannot solve for this problem, but why does it provide x-dependent eigenvalues?

EDIT: After bbgodfrey's answer, I'd like to go deep into this problem. The solutions are HeunT functions and the eigenvalues exist. In Dong 2019 you have an in-depth reference to this problem. Mathematica should be able to solve in version 12.1 since it has Heun triconfluent functions.

  • $\begingroup$ Thanks for the reference, the abstract of which is consistent with my answer below. Unfortunately, I do not have access to the article itself. Does it give symbolic solutions for the eigenvalues too, or only for the eigenfunctions? By the way, in my experience DSolve, on which DEigensystem depends, fails to solve many ODEs for which solutions are known.. What troubles me more, though, is that DEigensystem gives spurious answers for {x, -∞, ∞}, instead of returning unevaluated. Perhaps, it is time to report this as a bug. $\endgroup$
    – bbgodfrey
    Apr 15, 2020 at 12:45
  • $\begingroup$ @bbgodfrey They give symbolic solutions for the eigenfunctions, but the eigenvalues are computed numerically. I have the paper, I would not attach here but I can send it to you, if you're interested. I guess reporting a bug would be useful for other users. $\endgroup$
    – Fabio
    Apr 15, 2020 at 14:24
  • $\begingroup$ I am happy to report it as a bug and add the usual bug header, if you like. $\endgroup$
    – bbgodfrey
    Apr 15, 2020 at 15:00
  • 1
    $\begingroup$ ok, let's do it! :) $\endgroup$
    – Fabio
    Apr 15, 2020 at 15:18

2 Answers 2


Version 12.1 also give this strange result, which may be a bug.

DEigensystem[{-1/2 y''[x] + (-x^2/2 + x^4/4) y[x]}, y[x], {x, -∞, ∞}, 2]
(* {{Sqrt[-2 + x^2]/(2 Sqrt[2]), (3 Sqrt[-2 + x^2])/(2 Sqrt[2])}, 
   {E^(-((x^2 Sqrt[-2 + x^2])/(2 Sqrt[2]))), 
      2^(3/4) E^(-((x^2 Sqrt[-2 + x^2])/(2 Sqrt[2]))) x (-2 + x^2)^(1/4)}} *)

On the other hand, reducing the limits to

DEigensystem[{-1/2 y''[x] + (-x^2/2 + x^4/4) y[x]}, y[x], {x, -4, 4}, 2]

returns unevaluated. It is, however, easy to solve this problem numerically.

snn = NDEigensystem[{-1/2 y''[x] + (-x^2/2 + x^4/4) y[x]}, y[x], {x, -4, 4}, 4]; 
snn // First
Plot[Evaluate[snn // Last], {x, -4, 4}, PlotRange -> All]
(* {0.147275, 0.872551, 2.12949, 3.59638} *)

enter image description here Returning now to the symbolic solution problem posed in the question, consider

s = DSolveValue[-1/2 y''[x] + (-x^2/2 + x^4/4) y[x] == lamda y[x], y[x], x]
(* E^((x (3 - x^2))/(3 Sqrt[2])) C[1]^2 
       HeunT[-(1/2) - 2 lamda, -Sqrt[2], Sqrt[2], 0, -Sqrt[2], x] + 
   E^(-((x (3 - x^2))/(3 Sqrt[2]))) C[2] 
       HeunT[-(1/2) - 2 lamda, Sqrt[2], -Sqrt[2], 0, Sqrt[2], x] *)

The HeunT functions, newly defined in Version 12.0, are not bounded for large Abs[x], which may confuse Mathematica. However, if DEigensystem, cannot solve this problem, it should return unevaluated. As noted above, this may be a bug. The first eigenvalue can, however, be obtained by

sr1 = FindRoot[0 == D[s, x] /. {C[1]^2 -> 1, C[2] -> 1, x -> 4}, 
    {lamda, 0.15}, WorkingPrecision -> 45]
(* {lamda -> 0.147235140084093444055886856920046970641521561} *)

which agrees with the first eigenvalue determined by NDEigensystem above to four significant figures. (Note that FindRoot does not converge for smaller WorkingPrecision.) Likewise, the second eigenvalue can be obtained by

sr2 = FindRoot[0 == D[s, x] /. {C[1]^2 -> 1, C[2] -> -1, x -> 4}, 
    {lamda, 1}, WorkingPrecision -> 45]
(* {lamda -> 0.872261197867424491236128918505849507397539635} *)

Larger even and odd eigenvalues also are given by sr1 and sr2, respectively, with larger initial guesses for lamda. Substituting the values for {C[1]^2, C[2], lamda} into s then yields the same curves as above, up to a normalization factor. (Plot requires WorkingPrecison -> 30 for smooth curves.) So, with assistance Mathematica can obtain symbolic solutions for the eigenfunctions, although not for the eigenvalues.

I would welcome the readers' views on whether the behavior of DEigensystem is a bug. (Trace produces very lengthy but not particularly informative output.)

Addendum - Computation Using New Feature of DSolve 12.1

Documentation for DSolve 12.1 (under Scope) describes how to solve Sturm-Liouville problems. It can be applied to the present problem as follows.

newds = DSolveValue[{-1/2 y''[x] + (-x^2/2 + x^4/4) y[x] == lamda y[x], 
    y'[-4] == 0, y'[4] == 0}, y[x], x, Assumptions -> 0 < lamda < 4];

producing a lengthy Piecewise function that contains the eigenfunction with one constant of integration eliminated and a transcendental equation for lamda (not the actual eigenvalues in this case). This equation is extracted by newds[[1, 1, 2, 1]], which is solved without difficulty for the eigenvalues.

Table[FindRoot[newds[[1, 1, 2, 1]] /. C[1] -> 1, {lamda, n}, 
     WorkingPrecision -> 45], {n, .2, 3.2, 1}] // Values // Flatten // N[#, 6] &
(* {0.147235, 0.872261, 2.12798, 3.59109} *)

as expected.

  • $\begingroup$ Hi, I've added a reference to the analytical solutions of this problem. your analysis is interesting. Mathematica should solve this problem analytically and provide the correct eigenvalues, since the solutions are Heun triconfluent functions $\endgroup$
    – Fabio
    Apr 15, 2020 at 7:15
  • $\begingroup$ What is the reason to put c[1],c[2]$=\pm1$ in the code for sr1,sr2? I feel they can have different magnitude in general? $\endgroup$
    – xiaohuamao
    Oct 7, 2020 at 23:08
  • $\begingroup$ @xiaohuamao I chose C[1] = C[2] and C[1] = - C[2], because I knew that doing so would give symmetric and antisymmetric solutions, consistent with the numerical plot in my answer. To verify this, run FindRoot[{0 == D[s, x] /. x -> 4, 0 == D[s, x] /. x -> -4} /. C[1] -> 1, {{lamda, 0.15}, {C[2], 2}}, WorkingPrecision -> 45], which returns lambda as before and C[2} -> 1.. Similarly, FindRoot[{0 == D[s, x] /. x -> 4, 0 == D[s, x] /. x -> -4} /. C[1] -> 1, {{lamda, 1}, {C[2], 2}}, WorkingPrecision -> 45] returns the second eigenvalue as before and C[2} -> -1.. $\endgroup$
    – bbgodfrey
    Oct 8, 2020 at 1:13
  • $\begingroup$ Thanks. Now I see this is partially because the differential operator is parity even and the Sturm-Liouville theorem. $\endgroup$
    – xiaohuamao
    Oct 8, 2020 at 5:09
  • $\begingroup$ Dong 2019 is very nice paper, with now so many citations. So cannot its result that is integrated (apparently?) into DSolveValue[] be also integrated into DEgensystem? Also the other paper discusses how to improve Appendix A result from the paper [of Dong 2019] by using the Wronskian determinant: doi.org/10.1016/j.physleta.2020.126480 $\endgroup$ Jun 10, 2022 at 13:38

The question is still ill-posed. This is most probably not a sensation that the eigensystem is dependent on x. It is a sensation the n in the built-in is not chosen appropriately and as suggested by the Mathematica documentation. Prefer n=4 instead of n=2. The reason this is like the Laplacian a second-order differential equation with on x relying coefficient functions.

So the input is:

DEigensystem[-1/2 y''[x] + (-x^2/2 + x^4/4) y[x], 
 y[x], {x, -∞, ∞}, 4]


With my version 12.0.0 even n=2 will pose the results shown by others.

The main problem is that the interval between {-2,2} has no real values they are all imaginary. Therefore the system can not be properly solved on the complete interval entered in the DEigensystem.

So the result in the imaginary part {-2,2} is displayed with

Plot[Evaluate[Im@%28[[2]]], {x, -10, 10}, PlotRange -> Full]

graphics basis eigensystem

The real part:

Plot[Evaluate[%28[[2]]], {x, -10, 10}, PlotRange -> Full]

real part of the basis

This is in Mathematica 12.0.0. Increase the n option in yours might be appropriate.

  • 3
    $\begingroup$ To what does %28 refer? $\endgroup$
    – bbgodfrey
    Apr 15, 2020 at 21:21
  • 2
    $\begingroup$ dear @user2432923, honestly I don't understand your comment at all and how can be useful your suggestion in providing n=4 to find the correct solution. I also claim that the problem is well-posed. The eigenvalues cannot depend on the spatial variable, furthermore from a software perspective we have provided that symbol and the software should interpret this information accordingly. I think also that DEigensystem should give the correct result at least for the eigenfunctions since DSolve does. Likely, the eigenvalues cannot be computed in an analytical way... DEigensystem seems to be buggy. $\endgroup$
    – Fabio
    Apr 16, 2020 at 7:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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