Gauss-Seidel using Matlab advantages

Let's say that you have a grid of LxM points. On this grid you would like to solve, let's say, Laplacian(f) = 0. And let's say that you want to use the matlab abilities of working nicely with matrices and vectors. Solving it by Jacobi method is easy. How do you solve it by Gauss Seidel method?