  We present a simple method to solve spherical harmonics moment systems, such
as the the time-dependent $P_N$ and $SP_N$ equations, of radiative transfer.
The method, which works for arbitrary moment order $N$, makes use of the
specific coupling between the moments in the $P_N$ equations. This coupling
naturally induces staggered grids in space and time, which in turn give rise to
a canonical, second-order accurate finite difference scheme. While the scheme
does not possess TVD or realizability limiters, its simplicity allows for a
very efficient implementation in Matlab. We present several test cases, some of
which demonstrate that the code solves problems with ten million degrees of
freedom in space, angle, and time within a few seconds. The code for the
numerical scheme, called StaRMAP (Staggered grid Radiation Moment
Approximation), along with files for all presented test cases, can be
downloaded so that all results can be reproduced by the reader.
