clear all
close all
% -laplace(u) = f
% homogeneous Dirichlet b.c. x=-1,y=1
% homogeneous Neumann b.c.
% exact solution (x^2+2x-3)(y^2-2y-3)
dx = 1;
dy = 1;
m = [25,27];
x = linspace(-1,1,m(1));
y = linspace(-1,1,m(2));
h = 2./(m-1);
[X,Y] = meshgrid(x,y);
f = @(X,Y) -dx*2.*(Y.^2-2*Y-3)-dy*2.*(X.^2+2*X-3);
Ax = toeplitz([2/h(1)^2,-1/h(1)^2,zeros(1,m(1)-2)]);
Ay = toeplitz([2/h(2)^2,-1/h(2)^2,zeros(1,m(2)-2)]);
Ax(1,:) = [3/2/h(1)^2,-2/h(1)^2,1/2/h(1)^2,zeros(1,m(1)-3)];
Ay(m(2),:) = [zeros(1,m(2)-3),1/2/h(2)^2,-2/h(2)^2,3/2/h(2)^2];
Imx = speye(m(1));
Imx(1,:) = zeros(1,m(1));
Imy = speye(m(2));
Imy(m(2),:) = zeros(1,m(2));
A = kron(Imx,dy*Ay)+kron(dx*Ax,Imy);
b = f(X,Y);
Dindex = [1:m(2):m(2)*(m(1)-2)+1];
Dindex = [Dindex,(m(1)-1)*m(2)+1:m(1)*m(2)];
A(Dindex,:) = zeros(length(Dindex),prod(m));
A(Dindex,Dindex) = speye(length(Dindex))/(h(1)^2+h(2)^2);
b(Dindex) = 0;
% Neumann boundary conditions on x=-1 (y<1)
Nindex = [2:m(2)];%,2*m(2):m(2):(m(1)-1)*m(2)]
b(Nindex) = 0;
% Neumann boundary conditions on (-1,1)
% Neumann boundary conditions on y=1
A(m(2),:) = kron([1,zeros(1,m(1)-1)],...
		 [zeros(1,m(2)-3),1/2/h(2)^2,-2/h(2)^2,3/2/h(2)^2]);
Nindex = [m(2):m(2):(m(1)-1)*m(2)];
b(Nindex) = 0;
U = A\b(:);
figure
mesh(X,Y,reshape(U,m(2),m(1)))
xlabel('x')
ylabel('y')
errorInf = norm(U-reshape((X.^2+2*X-3).*(Y.^2-2*Y-3),prod(m),1),Inf)
