clear all
M = 30;
f = @(x) -5+10*(x(:,1)>0.5).*ones(length(x),1);
uD = @(x) x(:,1).*x(:,2);

I = eye(M);
e = ones(M,1);
T = spdiags([-e,4*e,-e],[-1,0,1],M,M);
A = sparse(M^2,M^2);

%*** create stiffness matrix
for j = 1:M
  A((j-1)*M+[1:M],(j-1)*M+[1:M]) = T;
end
for j = 2:M
  A((j-2)*M+[1:M],(j-1)*M+[1:M]) = -I;
  A((j-1)*M+[1:M],(j-2)*M+[1:M]) = -I;
end
%*** create rhs
t = linspace(0,1,M+2);
[xx,yy] = meshgrid(t(2:end-1),t(2:end-1));
b = 1/(M+1)^2*f([xx(:),yy(:)]);
b = reshape(b,M,M);
%*** incorporate inhomogeneous Dirichlet condition
[xx,yy] = meshgrid(t,t);
zz = zeros(size(xx));
zz( : , 1 ) = uD([xx( : , 1 ), yy( : , 1 )]);
zz( : ,end) = uD([xx( : ,end), yy( : ,end)]);
zz( 1 , : ) = uD([xx( 1 , : ); yy( 1 , : )]')';
zz(end, : ) = uD([xx(end, : ); yy(end, : )]')';

b( : , 1 ) = b( : , 1 ) + zz(2:end-1, 1 );
b( : ,end) = b( : ,end) + zz(2:end-1,end);
b( 1 , : ) = b( 1 , : ) + zz(  1 , 2:end-1);
b(end, : ) = b(end, : ) + zz(end , 2:end-1);
%*** compute solution
x = A\b(:);
zz(2:end-1,2:end-1) = reshape(x,M,M);

figure(1)
surf(xx,yy,zz)
figure(2)
spy(A)