Deriving the Schrödinger Equation from First Principles

why bother deriving Schrödinger’s equation?

If you’ve taken any quantum course, you’ve almost certainly seen this written on the board with very little ceremony:

i\hbar \frac{\partial \psi}{\partial t} = \hat{H}\,\psi

It’s presented as the equation of motion in quantum mechanics — but where does it actually come from?
Is it just “postulated”, or can we at least motivate it from ideas we already know: waves, energy, momentum?

Classically, every observable—position, momentum, energy—is a definite number.
Quantum mechanics replaces that picture: a particle is described by a wave function $\psi(\mathbf{r}, t)$
whose squared modulus gives a probability density. This post is about connecting those two worlds: starting from de Broglie’s wave hypothesis and simple plane waves seeing how the time-dependent Schrödinger equation almost forces itself on us.

who is this for?

This walkthrough is aimed at readers who:

Know basic classical mechanics and the energy–momentum relation for a non-relativistic particle,
Are comfortable with complex exponentials and derivatives of $e^{i(kx - \omega t)}$ ,
Have seen the Schrödinger equation written down, but never felt satisfied with “we just promote $p$ to an operator”.

We’re not doing a fully rigorous axiomatic derivation here.
Instead, the goal is to build physical and mathematical intuition for why the familiar form of the equation is so natural.

de Broglie waves: starting from a plane wave

de Broglie postulated that a free particle with momentum $p$ behaves like a plane wave:

\psi(x, t) = A \exp\!\left(i(kx - \omega t)\right), \qquad k = \frac{p}{\hbar},\quad \omega = \frac{E}{\hbar}

Taking partial derivatives:

\frac{\partial \psi}{\partial t} = -i\omega\,\psi = -\frac{iE}{\hbar}\,\psi \quad \Longrightarrow \quad E\,\psi = i\hbar\,\frac{\partial\psi}{\partial t}

\frac{\partial^2\psi}{\partial x^2} = -k^2\,\psi = -\frac{p^2}{\hbar^2}\,\psi \quad \Longrightarrow \quad p^2\,\psi = -\hbar^2\,\frac{\partial^2\psi}{\partial x^2}

These relations tell us how energy and momentum act on a plane wave: as time and space derivatives.

promoting classical energy to an operator

For a particle in a potential $V(x)$ , the classical energy is:

E = \frac{p^2}{2m} + V(x)

We promote $E$ and $p^2$ to operators acting on $\psi$ :

Classical quantity	Quantum operator
$E$	$i\hbar\,\partial_t$
$p$	$-i\hbar\,\partial_x$
$p^2$	$-\hbar^2\,\partial_x^2$

Substituting these into the classical expression $E = p^2/(2m) + V(x)$ gives:

i\hbar\,\frac{\partial\psi}{\partial t} = \left(-\frac{\hbar^2}{2m}\frac{\partial^2}{\partial x^2} + V(x)\right)\psi

The bracket on the right is the Hamiltonian operator $\hat{H}$ , giving us exactly the Schrödinger equation.

This argument is deliberately heuristic: we motivated the operator replacements using plane waves and then assumed they extend to more general wave packets and potentials. A more complete treatment comes from the axioms of quantum mechanics and symmetry arguments, but this route is often the most illuminating the first time you meet the equation.

a quick numerical check in Python

Let's verify the ground-state energy of the infinite square well, $E_n = \frac{n^2\pi^2\hbar^2}{2mL^2}$ ,
by finite-differencing the Hamiltonian:

python

import numpy as np
 
def infinite_square_well_energies(L: float, n_points: int = 500, n_states: int = 5):
    """
    Numerically diagonalise the 1-D particle-in-a-box Hamiltonian.
    Returns the first `n_states` eigenvalues in units of ℏ²/(2mL²).
    """
    hbar_2m = 1.0  # natural units
    dx = L / (n_points + 1)
    x = np.linspace(dx, L - dx, n_points)
 
    # Second-derivative finite-difference stencil
    diag = np.full(n_points,  2.0 / dx**2)
    off  = np.full(n_points - 1, -1.0 / dx**2)
    H = np.diag(diag) + np.diag(off, 1) + np.diag(off, -1)
    H *= hbar_2m
 
    eigenvalues = np.linalg.eigvalsh(H)
    return eigenvalues[:n_states]
 
L = 1.0
numerical = infinite_square_well_energies(L)
analytical = [(n**2 * np.pi**2) / (2 * L**2) for n in range(1, 6)]
 
for n, (num, ana) in enumerate(zip(numerical, analytical), start=1):
    print(f"n={n}: numerical={num:.4f}, analytical={ana:.4f}")

Running that gives sub-0.1 % error even at only 500 grid points you can see how the familiar $n^2$ spectrum of the infinite well emerges from a simple matrix eigenvalue problem.

normalisation and the Born rule

The wave function must be normalised:

\int_{-\infty}^{\infty} |\psi(x,t)|^2\,dx = 1