I'm currently doing a mathematica for physicists course and am struggling to solve a problem we were given!

I'm supposed to define an initial wavefunction in a harmonic oscillator (displaced, stretched or compressed ground state for example). Then I should propagate the wave function in time and use NDSolve to solve the resulting PDE. To keep the expressions simple, constants are set to 1.

Now I really don't know how to use mathematica well yet so I'm struggling to start this exercise. I first define an initial wavefunction.

ψinit[x_] := Exp[-x^2/2]

Then to propagate in time I can multiply it by the time evolution:

ψinit[x_, t_] := ψinit[x] Exp[-I*E_n*t/ħ]

But I do not understand how I'm supposed to get a PDE from this that I can solve.

Any help appreciated :) Thanks

  • $\begingroup$ You need to use the time-dependent Schrödinger equation (which is a PDE), not the time-independent one (which is an eigenvalue equation). What you wrote is the time propagation of an eigenstate of the time-independent Schrödinger equation, which will never get you any interesting dynamics. $\endgroup$
    – Roman
    Jan 26, 2019 at 21:37

1 Answer 1


A time-dependent wavefunction ψ[x,t] must satisfy the time-dependent Schrödinger equation $i\hbar\frac{\partial\psi(x,t)}{\partial t}=H\psi(x,t)$. For a harmonic oscillator with Hamiltonian $H=-\frac{\hbar^2}{2m}\frac{\partial^2}{\partial x^2}+\frac12m\omega^2x^2$ this would be

I*ħ*D[ψ[x,t],t] == -ħ^2/(2m)*D[ψ[x,t],{x,2}] + 1/2*m*ω^2*x^2*ψ[x,t]

Here's an example, using a finite calculation box and setting $\hbar=m=\omega=1$:

a = 5;
NDSolve[{I*D[ψ[x,t],t] == -1/2*D[ψ[x,t],{x,2}] + 1/2*x^2*ψ[x,t],
  ψ[x,0]==Exp[-x^2/2], ψ[-a,t]==ψ[a,t]==0}, ψ, {x,-a,a}, {t,0,10}]

(not being very careful with boundary conditions here).

To be more careful with the boundary conditions, let's use a beta distribution as a starting wavefunction: it is almost indistinguishable from a normal distribution, with the difference that it reaches zero precisely at the boundaries of the simulation box. For a Beta-Gaussian centered at $x=0$ with standard deviation $\sigma^2$ this would be

$$ \psi_0(x) \propto \left[1-\left(\frac{x}{a}\right)^2\right]^{\frac{1}{4} \left[\left(\frac{a}{\sigma}\right)^2-3\right]} $$

a = 5;
σ = 0.5;
ψ0[x_] = (1 - (x/a)^2)^(((a/σ)^2 - 3)/4);
NDSolve[{I*D[ψ[x, t], t] == -1/2*D[ψ[x, t], {x, 2}] + 1/2*x^2*ψ[x, t],
  ψ[x, 0] == ψ0[x], ψ[-a, t] == ψ[a, t] == 0}, ψ, {x, -a, a}, {t, 0, 10}]

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.