How many different ways can one solve an eigenvalue problem and plot its corresponding eigenfunctions in Mathematica? For example for Harmonic Oscillator? Which one is the most accurate one?

Thanks in advance.

  • $\begingroup$ You might look here: link $\endgroup$
    – eldo
    Commented Jun 13, 2014 at 9:48
  • $\begingroup$ Actually, I am not looking for Harmonic Oscillator. I am intersted how can one in general solve Eigenvalue Problem numerically by Mathematica. $\endgroup$
    – user14782
    Commented Jun 13, 2014 at 9:52
  • $\begingroup$ What sort of eigen problem are you trying to solve? Is it a partial differential equation or a set of ordinary differential equations? Some code from you would help. Have you looked at Eigensystems and its variants in Help? $\endgroup$
    – Hugh
    Commented Jun 13, 2014 at 10:20
  • $\begingroup$ it is an ordinary 1 dimensional second order differential equation. Imagine potential is 1/2x^2 and I want to obtain eigenvalues and plot eigenfunctions. I would like to learn the general procedure of solving such problems in Mathematica. $\endgroup$
    – user14782
    Commented Jun 13, 2014 at 10:23
  • $\begingroup$ @user14782 the method will depend on what Schroedinger equation you want to solve. Single particle or many-body? In the continuum or on a lattice? Is the Hamiltonian time-dependent or not? Even if you specify single-particle, in a continuum and time-independent, you can 1) discretise space and solve, 2) solve exactly part of the problem (eg the potential part), expand the rest in terms of that, and trunacate (1 is a special case of this) etc. Basically you need to be more specific, otherwise the answer to your question would be the bulk of modern computational physics. $\endgroup$
    – acl
    Commented Jun 13, 2014 at 10:39

2 Answers 2


Here is how you can solve the simple harmonic oscillator — i.e. quadratic potential — eigenvalue problem using Mathematica. For simplicity, I set all of the constants to unity.

Define the differential (Schrödinger) equation.

deqn = -(1/2) y''[x] + 1/2 x^2 y[x] == e y[x];

Solve the differential equation.

sol = DSolve[deqn, y, x][[1]]

  {y -> Function[{x}, 
    C[2] ParabolicCylinderD[1/2 (-1 - 2 e), I Sqrt[2] x] + 
    C[1] ParabolicCylinderD[1/2 (-1 + 2 e), Sqrt[2] x]]}

This is the general solution, parameterised by the eigenvalue e and two constants of integration C[1] and C[2].

We have to ensure that the solution is square integrable, so it had better behave itself as x goes to +Infinity and as x goes to -Infinity.

Do a series expansion about x = +Infinity.

y[x] /. sol // Series[#, {x, Infinity, 1}] & // Expand

(* E^(x^2/2) C[2] (expression1) + E^(-(x^2/2)) C[1] (expression2) *)

So, rather than sol, you must now use sol /. C[2] -> 0, so that C[2] = 0 suppresses the divergent E^(x^2/2) term.

Do a series expansion about x = -Infinity.

y[x] /. sol /. C[2] -> 0 // Series[#, {x, -Infinity, 1}] & // Expand

(* E^(x^2/2) (expression1a) + E^(-(x^2/2)) (expression2a) *)

Somehow you must arrange things so that expression1a suppresses the divergent E^(x^2/2) term. In this case expression1a contains the factor

suppress = -((I 2^(1/4 (1 - 2 e)) E^(-I e \[Pi]) Sqrt[\[Pi]])/ Gamma[1/2 (1 - 2 e)]);

If you plot suppress versus e you find that it has zeros at half-integer values of e — i.e. e = 1/2 + m where m = 0, 1, 2, ... — which gives you the familiar ladder of simple harmonic oscillator eigenvalues.

Plot[Abs[suppress], {e, 0, 5}]

(* graphics *)

So the (yet to be normalised) solution of the differential equation is finally

y[x] /. sol /. C[2] -> 0 /. e -> 1/2 + m // Simplify

(* C[1] ParabolicCylinderD[m, Sqrt[2] x] *)

Verify that this is indeed a solution of the differential equation.

deqn /. sol /. C[2] -> 0 /. e -> 1/2 + m // FullSimplify[#, Assumptions -> m \[Element] Integers] &

(* True *)

Variants of this approach can be used solve for the eigenfunctions and eigenvalues of other (sufficiently simple) potentials.

  • $\begingroup$ It's Great! Thanks! $\endgroup$
    – user14782
    Commented Jun 13, 2014 at 17:12
  • $\begingroup$ @user14782 - you should click the up- and down-buttons of Stephen Luttrel's answer. $\endgroup$
    – eldo
    Commented Jun 13, 2014 at 18:00
  • $\begingroup$ You are right! I was not familiar with that. $\endgroup$
    – user14782
    Commented Jun 13, 2014 at 18:29
  • $\begingroup$ @user14782 ideally you only click the up button if you find the answer helpful, of course :) $\endgroup$
    – acl
    Commented Jun 13, 2014 at 19:42

You can solve the characteristic polynomial of your matrix (equation(s)):

CharacteristicPolynomial[ m, \[Lambda]]
Solve[0 == %, \[Lambda]]

or just use this function

  • $\begingroup$ I would like to solve Schrodinger differential equation and plot the eigenfunctions. $\endgroup$
    – user14782
    Commented Jun 13, 2014 at 10:15