bsize = 1.;
np = 20;
niter = 10000;

tic
[u, xx, yy] =  myfd2dpoisson(@rhs, @bound, bsize, np, niter);
toc

[x, y] = meshgrid(xx, yy);  % mesh for plotting

figure(1);
title('Numerical solution');
surf(x, y, u);

uexact = zeros(np+2, np+2);
for i = 1:np+2
    x1 = xx(i);    
    for j = 1:np+2
        y1 = y(j);
        uexact(i, j) = bound(x1, y1);
    end
end
   
figure(2);
title('Exact solution');
surf(x, y, uexact);

figure(3)
title('Error');
surf(x, y, u-uexact);

norm(u-uexact, inf)

function f = rhs(x, y)
    % x and y expected to be are scalars
	f = (x^2 + y^2)*exp(x*y);
end

function g = bound(x, y)
    % x and y expected to be are scalars
 	g = exp(x*y);
end
