Simplest way is to use Gauss-Seidel and solve it iteratively. Laplace's equation is particularly easy here since (for an equally-spaced grid) each nodal value is just the average of the four nodes around it. Can usually be accelerated with over-relaxation (PSOR).
Faster (but more complex) iterative method is ADI (alternating-direction implicit), with repeated sweeps of a tri-diagonal solver.