function [A,b] = lap2dd

L = 1.0;
W = 1.0;
l =10;
m = 10;
K = 2;
C = 0.01;
c = 0.01;
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 i=1:nx
    A(i,:) = 0;
    A(i,i) = 1;
    b(i) = Ubc((i-1)*dx,0);
    
    A((ny-1)*nx+i,:) = 0;
    A((ny-1)*nx+i,(ny-1)*nx+i) = 1;
    b((ny-1)*nx+i) = Ubc((i-1)*dx,1);
end

for j=1:ny
    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);
    
    A(j*nx,:) = 0;
    A(j*nx,j*nx) = 1;
    b(j*nx) = Ubc(1,(j-1)*dy);
end



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