replace conv numerical diffusion with matrix mul

this allows neumann BCs to be built into the operation - chose 2nd order
differences with derivative η=0

a future improvement would be to choose η to be +/- 1

fixes #2