function v = Ax(u);
% on input: u is a N^2 vector

n = length(u);
N = sqrt(n);
U = reshape(u, [N N]);

V = 4.0 * U;
V( :     ,2:end  ) = V( :     ,2:end  ) - U( :     ,1:end-1);
V(2:end  , :     ) = V(2:end  , :     ) - U(1:end-1, :     );
V(1:end-1, :     ) = V(1:end-1, :     ) - U(2:end  , :     );
V( :     ,1:end-1) = V( :     ,1:end-1) - U( :     ,2:end  );

v = V(:); % = reshape(V, [n 1]);

