%% Checking a solution that should be exact

f = @(x,y) 6*y-2;
g = @(x,y) 1+x.^2-y.^3;
M=10;
[U,X,Y]=fdPoisson(f,g,M);
exact = g(X,Y);
disp('Error: '), disp(max(abs(exact(:)-U(:))));
surf(X,Y,U)

%% Checking if order is two with two grids

f = @(x,y) -2*cos(x.^2+y)+(4*x.^2+1).*sin(x.^2+y);
g = @(x,y) sin(x.^2+y);

M=10;
[U,X,Y]=fdPoisson(f,g,M);
exact = g(X,Y);
disp('Error: '), disp(max(abs(exact(:)-U(:))));

M=20;
[U,X,Y]=fdPoisson(f,g,M);
exact = g(X,Y);
disp('Error: '), disp(max(abs(exact(:)-U(:))));

%% Checking if order is two with many grids

f = @(x,y) -2*cos(x.^2+y)+(4*x.^2+1).*sin(x.^2+y);
g = @(x,y) sin(x.^2+y);

listM=2:2:20;
%listM=2.^(2:8);
h=1./listM;
error=[];
for M=listM
    [U,X,Y]=fdPoisson(f,g,M);
    exact = g(X,Y);
    error=[error max(abs(exact(:)-U(:)))];
end
% Log-log plot (experiments go from right-to-left)
subplot(121)
loglog(h,error,'-o'), hold on
loglog(h,h.^2,'k')   % a line of the expected slope
% By-hand log-log
subplot(122)
plot(-log10(h),-log10(error),'-o'), xlabel('-log_{10}(h)'), ylabel('-log_{10}(error)')