(see e.g., Numerical Recipes in C , Ch 19.5)
We want to solve the Laplace equation
$$ \boldsymbol{\nabla}^2 u = 0 $$with given boundary conditions.
Rewrite as a diffusion equation (with $D=1$):
$$ \boldsymbol{\nabla}^2 u = \frac{\partial u}{\partial t} $$The equilibrium solution of the diffusion equation (i.e., after very long $t$ when relaxed to equilibrium)
$$ \lim_{t\rightarrow+\infty} u(\mathbf{x}, t) \equiv u(\mathbf{x}) $$is the solution of the boundary value problem.
For a visualization of how the solution for the 2D "wire on the box" problem relaxes towards the static solution of the Laplace equation, see the SOR movie: the process looks like the 1D heat diffusion problem solution, with the y-axis of the electrostatic problem playing the role of time in the 1D diffusion problem.
Why do these two problems look similar?
Set up the finite difference scheme for the 2D diffusion equation, using $u(x, y, t) = u(i\Delta, j\Delta, n\Delta t) = u_{ij}^n$:
$$ \begin{split} \frac{1}{\Delta^2}\left[(u_{i+1,j}^n -2u_{ij}^n + u_{i-1,j}^n) + (u_{i,j+1}^n -2u_{ij}^n + u_{i,j-1}^n)\right] \\ = \frac{1}{\Delta t}(u_{ij}^{n+1} - u_{ij}^n) \end{split} $$In 2D, this leap frog algorithm is stable for
$$ \frac{\Delta t}{\Delta ^2} \le \frac{1}{4} $$Let's assume that we can push our solution to the edge of stability and use $\frac{\Delta t}{\Delta ^2} = \frac{1}{4}$. Then the finite difference scheme becomes
$$ \frac{1}{4} \left(u_{i+1,j}^n + u_{i-1,j}^n) + (u_{i,j+1}^n + u_{i,j-1}^n\right) - u_{ij}^n = u_{ij}^{n+1} - u_{ij}^n $$or when collecting terms, the Jacobi algorithm re-emerges
$$ u_{ij}^{n+1} = \frac{1}{4} \left(u_{i+1,j}^n + u_{i-1,j}^n + u_{i,j+1}^n + u_{i,j-1}^n\right) $$where you now interpret $n$ as the step number in the iterations towards convergence.
(For real fields $u(\mathbf{x}, t)$.)