function [A,b] = lap2dn

L = 1.0;
W = 1.0;
l =10;
m = 10;
K = 2;
C = 0.0;
c = 0.0;
Usur = 0.0;

dx = L/l;
dy = W/m;

nx = l+1;
ny = m+1;
npts = (nx)*(ny);

A = spalloc(npts,npts,5*npts);

D = ones(npts,1)*(2*K*(1/(dx^2) + 1/(dy^2)) + C);
Dx = ones(npts,1)*K/(dx^2);
Dy = ones(npts,1)*K/(dy^2);

A = A + spdiags([-Dy -Dx D -Dx -Dy],[-nx -1 0 1 nx],npts,npts);

b = ones(npts,1)*C*Usur;

for j=1:ny-1
    A(j*nx,j*nx+1) = 0;
end
for j=2:ny
    A((j-1)*nx+1,(j-1)*nx) = 0;
end

for i=1:nx
    A(i,i) = A(i,i) - 2*c/dy;
    A(i,i+nx) = A(i,i+nx) - K/(dy^2);
    b(i) = b(i) - 2*c*Usur/dy;
    A((ny-1)*nx+i,(ny-1)*nx+i) = A((ny-1)*nx+i,(ny-1)*nx+i) + 2*c/dy;
    A((ny-1)*nx+i,(ny-1)*nx+i-nx) = A((ny-1)*nx+i,(ny-1)*nx+i-nx) - K/(dy^2);
    b((ny-1)*nx+i) = b((ny-1)*nx+i) + 2*c*Usur/dy;
end

for j=1:ny
    A(j*nx,j*nx) = A(j*nx,j*nx) + 2*c/dx;
    A(j*nx,j*nx-1) = A(j*nx,j*nx-1) - K/(dx^2);
    b(j*nx) = b(j*nx) + 2*c*Usur/dx;
    A((j-1)*nx+1,:) = 0;
    A((j-1)*nx+1,(j-1)*nx+1) = 1;
    b((j-1)*nx+1) = Ubc(0,(j-1)*dy);
end



function Ubc = Ubc(x,y)
h = 10;
if abs(y - 0.5) > 0.1 - 1e-8
    Ubc = 0;
else
    Ubc = h*(0.1 - abs(y - 0.5));
end



% function Ubc = Ubc(x,y)
% h = 10;
% Ubc = h*(0.5 - abs(y - 0.5));
