You create matrices that represent the operator. For instance, suppose you had a 1D function and were interested in the evenly spaced points (x1, x2, x3, x4, x5), with the function values at these points being

  y1
  y2
  y3
  y4
  y5

If the differential equation had a first derivative in it, you'd construct something like this, multiplied by some constant (e.g. 1/(2*dx) for an unscaled derivative):

  ? ? ? ? ?
  -1 0 1 0 0
  0 -1 0 1 0
  0 0 -1 0 1
  ? ? ? ? ?

So the derivative at each element is defined by the difference of the next and previous elements. Multiplying the column of function values by this gives you the derivatives. This doesn't work for the first and last element, and in fact you'll usually modify these rows depending on what boundary condition is needed, so I've just left them filled in with "?".

For a solver, you don't know what the y values actually are, so you construct a column that corresponds to the right side of the differential equation. For instance, if the equation was something like the trivial dy/dx = c and you were using the operator above the column would be

  ?
  c
  c
  c
  ?

with the first and last values to be filled in based on the boundary conditions. You then left matrix divide that by the operator matrix (i.e. multiply it by the inverse of the operator matrix). That gives the solution to the equation.

This is just a simple example and in practice the matrices will be larger and more complex.

x is f in the article

A is the five-point laplacian (a symmetric matrix with five non-zero diagonals), as implemented in the "step" function in the article, let's call it L

b is the datum h

If you set-up the sparse matrix L and the vector b (with the same code that performs the iterations in the article), then you can solve Poisson equation "L*f=h" by doing "f=L\h".