I have a problem with unexpected results of the matrix*vector multiplication fvMatrix<scalar>::Amul().

The code below is supposed to calculate the laplacian of the field psi, where psi=y*y is the square of the y-coordinate.
Because d/dy d/dy (y*y) == 2 , I expect to calculate result==2, but result varies between -1e-5 and 5e-7 on a 20x20x1 mesh. The result is uniform except for one layer next to a wall.

I think there is a serious flaw about how I call the Amul() function. Could you please check for an obvious mistake?
(The code can be run for example on the mesh from the icoFoam/cavity tutorial. I have defined fixedValue boundary conditions for result.
I would like to solve this issue, because it is needed to verify parts of my CFD code.)