I've just finished learning the physics behind the problem and would like to write a program in C++ than can solve the problem. I'm actually stuck at the start. I've quite a bit of research, the problem is there is not too many examples of code used to solve the problem.

I'm going to solve the problem using finite-difference form.

$$\dfrac{d^2\psi}{dx^2}\approx \dfrac{\psi_{n+1}+\psi_{n-1}-2\psi_{n}}{(\Delta x)^2}$$

Which allows us to rearrange in the form

$$\psi_{n+1}=2\psi_{n} - \psi_{n-1}-2(\Delta x)^2(E-V_n)\psi_n$$

Using the even-parity solution, we have

$$\psi(0)=1 \quad \quad \psi'(0)=0$$

at $n=0 \implies x=0$

Letting $m=1$ and $\hbar = 1$ I know that at the ground state $E_g =\frac{\pi^2}{8}$

I'm struggling to write code that can help me find $\psi_n$ for values of $n$.

Can I assume that $V_0 = 0$ since we want the wave function to be inside the well?

I'm not asking for someone to write the code for me, I'd just like tips on what I need to define and which method I should use.

Thanks in advance

EDIT: This is what I have so far

#include <iostream>
#include <string>
#include <stdio.h>
#include <unistd.h>
#include <math.h>
#include <stdlib.h>
#include <stdarg.h>
#include <assert.h>

using namespace std;

const double PI = 3.14159265358979323846264338327950;
int main() {

    cout < "1D Particle-in-a-Box"\n;

    double psi0, dpsi0, N, dx, x_end, E;
    int number_steps, nr;

    double value_x [number_steps];
    double psi [number_steps];
    double V [number_steps];

    //read parameters (even-parity)
    cout < "psi0 = ";
    cin > psi0;
    cout < "dpsi0 = ";
    cin > dpsi0; 

    dx = x_end / number_steps;

    value_x [0] = 0;
    psi [0] = psi0;
    V [0] = 0.0
    //finite-difference form

    while (nr < number_steps)
        value_x [nr+1] = dx*(nr+1); 
        psi[nr+1] = 2*psi[nr] - psi[nr-1] - (2*(dt*dt))*(E-V[nr])*(psi[nr]);
        nr = nr + 1;
  • $\begingroup$ What method do you want to use? Check this answer, for example $\endgroup$
    – nicoguaro
    Commented Nov 21, 2017 at 21:19
  • $\begingroup$ Shooting method and runge kutta since that is the easiest apparently. $\endgroup$ Commented Nov 22, 2017 at 18:25
  • $\begingroup$ Finite difference is the easiest, in my opinion. $\endgroup$
    – nicoguaro
    Commented Nov 23, 2017 at 2:31
  • 1
    $\begingroup$ I don't understand your boundary conditions. Is the box bounded by infinite barriers? If so, assuming the box is between $-L$ and $L$ , you must have on the one hand $\psi(-L)=\psi(L)=0$. On the other hand, the depth of the box being arbitrary, you can put $V(x)=0\ ;\ x \in \{-L;L\}$ no matter what. The solution being even-parity tells you that it is symmetrical with regards to the center of the box, not the value at any $x$. This depends on the potential. Plus, in the ground state your solution is a wave. Saying $n=0 \rightarrow x=0$ makes no sense in this problem. $\endgroup$
    – G.Clavier
    Commented Dec 29, 2017 at 11:56
  • $\begingroup$ @Patrick Moloney I think you should focus more on the physics of the basic problem before moving to try to program it. As pointed out by G. Clavier, It seems you have some misunderstandings of the underlying problem. I would also recommend looking over this paper as one way of solving generic 1d potential problems:google.com/url?sa=t&source=web&rct=j&url=http://… $\endgroup$
    – Tyberius
    Commented Feb 26, 2018 at 3:20

2 Answers 2


I'm not too familiar with QM, but I would use FDM. You could treat it as a steady-diffusion equation with a source term. You would have

$$ -\frac{h^{2}}{2 m} \frac{\psi_{i-1} - 2 \psi_{i} + \psi_{i+1}}{\Delta x^{2}} + V \psi_{i} = E \psi_{i} $$

You can rearrange this to be

$$ -\frac{h^{2}}{2 m} \frac{\psi_{i-1} - 2 \psi_{i} + \psi_{i+1}}{\Delta x^{2}} = (E-V)\psi_{i} $$

You can use TDMA (Tridiagonal Matrix Algorithm) to solve this ODE as the left-hand side is a tridiagonal matrix and the right-hand side is a forcing function. Therefore for the tridiagonal matix you have

$$ a_{i} = -\frac{h^{2}}{2m \Delta x^{2}}\\ b_{i} = \frac{h^{2}}{m \Delta x^{2}}\\ c_{i} = -\frac{h^{2}}{2m \Delta x^{2}}\\ d_{i} = (E - V) \psi_{i} $$

For the boundary condition $\psi(0) = 0$ you make the adjust to the matrix (assuming 1 is your starting index)

$$ a_{1} = 0\\ b_{1} = 1\\ d_{1} = 0 $$

For the Neumann BC $\psi'(0) = 0$ you will need to evaluate $\psi''$ at $i=1$ and then use a first order approximate $\psi'_{\partial D}$ and then plug in your first order approximate to $\psi''$ evaluate at the boundary.

EDIT: I noticed you said you wanted to use RK4, but I'm confused as you would do that for a Time Dependent Schrodinger Equation TDSE as opposed to TISE. If you are solving TDSE then you could use RK4 for $\frac{\partial \psi}{\partial t},$ however since type of PDE is parabolic I would suggest using implicit time stepping, such as Crank-Nicholson. In the transient case you would have

$$ -\frac{h^{2}}{2 m} \frac{\psi_{i-1} - 2 \psi_{i} + \psi_{i+1}}{\Delta x^{2}} + V_{i} \psi_{i} = i h \frac{\psi^{n+1}_{i} - \psi^{n}}{\Delta t} $$

This would give you a $O(\Delta x^{2}, \Delta t)$ approximation.

In explicit form this would give you

  • $\begingroup$ The equation is similar to diffusion considering imaginary time $\endgroup$
    – nicoguaro
    Commented Nov 25, 2017 at 4:41
  • $\begingroup$ I agree, to me it seems like it is a parabolic PDE with a forcing function. $\endgroup$
    – Simon
    Commented Nov 27, 2017 at 19:05

Easiest solution: set up a grid $x_0, \ldots, x_N$ (with spacing $\Delta x$) large enough that your wavefunction (approximately) vanishes at the boundary. Then set up the (symmetric tridiagonal) Hamilton matrix

$$H_{ij}= \begin{cases} V(x_i) + \frac{1}{\Delta x^2} && i=j\,,\\ -\frac{1}{2\Delta x^2} && |i-j| = 1\,, \\ 0 && \text{otherwise,} \end{cases} $$

where $V(x)$ is the potential and $\left\{-\frac{1}{2\Delta x^2}, \frac{1}{\Delta x^2}, -\frac{1}{2\Delta x^2}\right\}$ is the usual three-point finite-difference approximation to the kinetic energy operator $-\frac12 \frac {\partial^2}{\partial x^2}$.

Diagonalize this matrix (e.g. using Eigen) to obtain $M = U E U^T$. The $i$-th column of $U$ contains the eigenstate of you system, and the eigenvalue $E_i$ is the corresponding eigenenergy.

That's as simple as it can get but at the same time very accurate for most one-dimensional problems.

From your question it seems you want to use the shooting method or similar things. Don't do that. There is no reason in applying the shooting method to a linear one-dimensional problem: it's inefficient and not very stable.

Here is some code to try yourself:

#include <iostream>
#include <Eigen/Dense>

int main()
    int N = 200;
    double x0 = -10.0;
    double xN = 10.0;     
    auto potential = [] (double x) {return 0.5*x*x;};

    //set up the Hamiltonian
    Eigen::MatrixXd H = Eigen::MatrixXd::Zero(N,N);
    double dx = (xN - x0) / (N-1);

    for(int i=0; i<N; ++i)
        double x= x0 + i*dx;

        if(i>0) { H(i-1,i) = -0.5/(dx*dx);  }
        H(i,i) = potential(x) + 1.0/(dx*dx);
        if(i<(N-1)) { H(i+1,i) = -0.5/(dx*dx);  }

    //... and diagonalize it
    Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eigensolver(H);
    std::cout<< eigensolver.eigenvectors().transpose() << std::endl;
  • $\begingroup$ Just to clarify, is the expression for the kinetic energy obtained by using the 3 point finite difference approximation for a 2nd derivative? $\endgroup$
    – Tyberius
    Commented Mar 4, 2018 at 22:33
  • $\begingroup$ @Tyberius: yes. $\endgroup$
    – davidhigh
    Commented Mar 5, 2018 at 8:00

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.