// Let N(t,i,j) represent the neutron density at time t at grid position (i,j).
// Initial conditions
for all grid positions (i,j)
N(0,i,j) = 0 if not at center of grid
N(0,i,j) = 10^7 at center of grid
for each time step t,
for all grid positions (i,j) (except the outside boundary)
N(t+1,i,j) = stencil(i,j)
// Boundary conditions
N(t+1,boundary) = 0
Inspect the simulation if desired.
endfor