
% Numerical approximation to the solution of Poisson's equation 
%
%   U_xx + U_yy = pf(x,y)
%
% over the square [0,a]x[0,a] with Dirichlet boundary conditions  
%
%   U(x,y) = gf(x,y)
%
% Uses a uniform mesh with (n+2)x(n+2) total points 
% (i.e, n x n interior grid points).

a = 1;
n = 39;

pf = @(x,y) -(x.^2 + y.^2).*sin(x.*y);
gf = @(x,y) sin(x.*y);

tic();
[u, x, y] = myfd2poisson(pf, gf, a, n);
toc()

% Plot solution
figure;
surf(x, y, u);
h = x(1,2) - x(1,1);
title(strcat('Numerical Solution to Poisson Equation, h=', num2str(h)));
xlabel('x');
ylabel('y');
zlabel('u(x,y)');

% Plot error
figure;
surf(x, y, u - gf(x,y));
h = x(1,2) - x(1,1);
title(strcat('Numerical Solution to Poisson Equation, h=', num2str(h)));
xlabel('x');
ylabel('y');
zlabel('error');

% Systematic error plot

np = 16;
h = zeros(np, 1);
errs = zeros(np, 1);

for i = 1:np
    n = i*8;
    [u, x, y] = myfd2poisson(pf, gf, a, n);
    h(i) = x(1,2) - x(1,1);
    errs(i) = norm(u - gf(x,y), inf);
end
figure;
loglog(h, errs, 'r.-');
hold on
loglog(h, h/(np*100), '--');
grid on
title('Solution error vs discretization size');
xlabel('h');
ylabel('error');
legend('numerical error', 'h^1 guideline')
