%% UnitSqFD
%  
% solves a simple poisson equation 
%
% $\Delta u = -1$ in the unit square,
% $u = 0$ on boundary

%% Parameters
n = 10;
x = linspace(0,1,n+2);
y = x;
[X,Y]=meshgrid(x,y);  % necessary for graphics

%% Exact Solution
% by the method of Greens functions, in form of a series
Z = zeros(size(X));
for m=0:15
    zm = 2*m+1;
    s = sin(zm*pi*X)/zm^3;
    ch = 1-cosh(zm*pi*(Y-0.5))./cosh(zm*pi*0.5);
    Z = Z + s.*ch;
end
Z = Z*4/pi^3;
surf(X,Y,Z)

% %% Finite Diferences
% % using the built-in command for the discrete Laplacian
% h = 1/(n+1);
% A = gallery('poisson',n);
% b = h^2*ones(n*n,1);
% x = A\b;
% Zinner = zeros(n);
% Zinner(:) = x(:);
% Zh = zeros(size(Z));
% Zh(2:end-1,2:end-1) = Zinner;
% surf(X,Y,Z)

% %% Solution Error 
% surf(X,Y,Zh-Z)





