function ShowPoisson()
clc
close all
% Set up the problem...
m = 60;
n = 80;
A = TheLaplacian(m,n);
h = .01;
b = zeros(m*n,1);
b(1:n) = fTop(n,h);
b(m*n-n+1:m*n) = fBot(n,h);
b = b*h^2;
% Solution via direct solve...
[G,err,P] = chol(A,'lower');
uStar = P'*(G'\(G\(P'*b)));
uStar = A\b;
ShowA(reshape(uStar,n,m)')
title('Solution Via Cholesky','fontsize',16)
pause(5)

% The Players for Jacobi and Gauss-Seidel...
L = tril(A,-1); D = diag(diag(A)); U = L';

% Jacobi...
MJ = D; NJ = -(L+U);
MG = D+L; NG = -U;

nSteps = 1000;
% nSteps with starting vector = 0...
u = zeros(m*n,1);
for k=1:nSteps
u = MJ\(b + NJ*u);
if rem(k,5)==0
ShowA(reshape(u,n,m)')
home
disp('Jacobi')
disp(sprintf(' Iterations = %1d',k))
disp(sprintf(' Error = %7.4f',norm(uStar-u)))
end
end
title('Solution Jacobi','fontsize',16)
figure
ShowA(reshape(uStar,n,m)')
title('Solution Via Cholesky','fontsize',16)
pause
close all
% Gauss-Seidel
% nSteps with starting vector = 0...
u = zeros(m*n,1);
for k=1:nSteps
u = MG\(b + NG*u);
if rem(k,5)==0
ShowA(reshape(u,n,m)')
home
disp('Gauss-Seidel')
disp(sprintf(' Iterations = %1d',k))
disp(sprintf(' Error = %7.4f',norm(uStar-u)))
end
end
title('Solution Gauss-Seidel','fontsize',16)
figure
ShowA(reshape(uStar,n,m)')
title('Solution Via Cholesky','fontsize',16)

function z = fTop(n,h)
x = (1:n)'*h;
c = n*h/4;
z = 5*exp(-10*(x-c).^2);

function z = fBot(n,h)
x = (1:n)'*h;
c = n*h/1.3;
z = 5*exp(-20*(x-c).^2);

function ShowA(A)
% Graphically displays the matrix A using the 'hot' colormap
colormap('hot')
pcolor(A)
shading interp
caxis([min(min(A)) max(max(A))])
axis ij equal off
shg