% Solves the Poisson equation with Dirichlet boundary conditions in a unit 
% square given by 
%
%   u_{xx} + u_{yy} = \exp xy        (x,y) \in [0,1]^2,
%            u(0,y) = 5 y (y - 1)       y  \in [0,1],            
%            u(1,y) = - y ((y-1)^4)     y  \in [0,1],          
%            u(x,0) = 0.5 sin 6\pi x  x    \in [0,1],            
%            u(x,1) = sin 2\pi x      x    \in [0,1].
%
% Compares Gaussian eleimination (Matlab's backslash division)and fast sine 
% transfom methods.

clear, clf

x = linspace(0,1,128+1);  Dx = min(diff(x));
y = linspace(0,1,64+1)';  Dy = min(diff(y));
[X Y] = meshgrid(x,y);
f = exp(X.*Y);
f(:,1) = 5.*y.*(y-1);        % east
f(:,end) = -y.*((y-1).^4);   % west
f(1,:) = .5.*sin(6*pi.*x);   % south
f(end,:) = sin(pi*2.*x);     % north

subplot(221)
tic, u = ps2d(f,Dx,Dy); et = toc;
mesh(X,Y,u)
xlabel('x'), ylabel('y'), zlabel('u')
title({'Gaussian elimination';['(elapsed time = ' num2str(et) ')']})

subplot(222)
tic, u = fps2d(f,Dx,Dy); et = toc;
mesh(X,Y,u)
xlabel('x'), ylabel('y'), zlabel('u')
title({'fast sine transform';['(elapsed time = ' num2str(et) ')']})