It turns out that numerical solutions to Laplace's equation (del squared phi equals zero) are very easy to code.
phii,j = (phii+1,j+phii-1,j+phii,j+1+phii,j-1)/4
Repeat until phi doesn't change much. I've played with this algorithm before just for fun, but I didn't realize it could be applied to solve real problems.
This is a specific case of a general method for solving a linear system, where you just iteratively satisfy each row of the system -- if you overwrite values as you go, it's called the Gauss Seidel method; otherwise, it's called the Jacobi method.
If you're willing to generalize a bit, the idea shows up all over the place -- in resolving collisions for a simulation, you'll often see people go to each collision and locally resolve it, iterate a few times, and hope the system converges to a globally valid state. And for doing inverse kinematics with the simple "Cyclic Coordinate Descent" method, you go to each joint in order, solve for the locally best angle at that joint assuming none of the others move, and repeat until convergence.
It's usually not the most efficient solution, but it's easy to code and understand and often works surprisingly well