function [A] =diff2d(n,K,a,h)
% N x N  sparse Laplacian
%n = 3;
ih2 = 1/(h^2);
nzmax = 5*n^2; % assume 5 entries per row, i.e. ignores boundary conditions
rowind = [1:nzmax];
colind = rowind;
vals = rowind;
k= 1;
row = 1;

    for iy = 1:n
        yoff = (iy-1)*n;
        for ix = 1:n
            indx = yoff +ix;
            if(iy > 1)
                rowind(k) = row;   val(k) = -K*ih2;
                colind(k) = indx - n; k = k+1;
            end
            if(ix > 1)
                rowind(k) = row;   val(k) = -K*ih2;
                colind(k) = indx - 1; k = k+1;
            end
            rowind(k) = row;   val(k) = 4*K*ih2+a;
            colind(k) = indx;  k = k+1;
            if(ix < n)
                rowind(k) = row;   val(k) = -K*ih2;
                colind(k) = indx + 1; k = k+1;
            end
            if(iy < n)
                rowind(k) = row;   val(k) = -K*ih2;
                colind(k) = indx + n; k = k+1;
            end
            row = row +1;
        end
    end
nzmax = k-1;
rows = rowind(1:nzmax);
cols = colind(1:nzmax);
vals = val(1:nzmax);
A = sparse(rows,cols,vals,n^2,n^2,nzmax);