clear all
close all
% -laplace(u) = f
% homogeneous Dirichlet b.c.
% exact solution (x^3-x)*(y^2-1)
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*6*X.*(Y.^2-1)-dy*2*(X.^3-X);
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)]);
A = kron(speye(m(1)),dy*Ay)+kron(dx*Ax,speye(m(2)));
b = f(X,Y);
Dindex = [1:m(2)];
for i = 2:m(1)-1
  Dindex = [Dindex,m(2)+(i-2)*m(2)+1,(i-2)*m(2)+2*m(2)];
end
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);
%A(Dindex,Dindex) = 1e18*speye(length(Dindex))/(h(1)^2+h(2)^2);
b(Dindex) = 0/(h(1)^2+h(2)^2);
%b(Dindex) = 1e18/(h(1)^2+h(2)^2);
U = A\b(:);
mesh(X,Y,reshape(U,m(2),m(1)))
xlabel('x')
ylabel('y')
figure
mesh(X,Y,(X.^3-X).*(Y.^2-1))
xlabel('x')
ylabel('y')
errorInf = norm(U-reshape((X.^3-X).*(Y.^2-1),prod(m),1),Inf)
