%% BVP / Elliptic PDE examples (Steady-State Heat Eqn.)
% Solving Laplace equation Txx+Tyy=0 subject to boundary conditions (9.1)

%Distance in Centimeters
w=10;       %Width
h=15;       %Height
T0=100;     %Temperature Peak
dx=2.5;    
dy=dx;
x=[0:dx:w]; 
y=[0:dy:h];
[X,Y]=meshgrid(x,y);
T=T0.*sinh(pi*Y./w).*sin(pi*X./w)/sinh(pi*h/w)

%Plot Solution
figure(1);
imagesc(x,y,T); 
axis xy;
xlabel('Width (cm)');
ylabel('Height (cm');
title('Analytical Solution: Temperature (C)');
colorbar;

%Apply 2D PDE solution, directly
Td=zeros(size(T));
m=length(y)-2;
n=length(x)-2;
o=m*n;%(length(y)).*(length(x));
C=zeros(1,o)';
diagonal=sparse(1:o,1:o,-4*ones(1,o),o,o);
suprdiag=sparse(2:o,1:o-1,ones(1,o-1),o,o);
for i=1:o
    C(i)=-T0*(mod(i,m)==0).*(sin(pi*i*dx/(w*m)));
    if ((mod(i,m)==0) && (i<o))
        suprdiag(i+1,i)=0;
    end;
end;
banddiag=sparse(m+1:o,1:o-m,ones(1,o-m),o,o);
smatrix=suprdiag+diagonal+suprdiag'+banddiag+banddiag';
fmatrix=full(smatrix);
Tdd=fmatrix\C;
Td(2:end-1,2:end-1)=reshape(Tdd,m,n)
Td(end,:)=T(end,:);

%Plot Results
figure(2);
imagesc(x,y,Td); 
axis xy;
xlabel('Width (cm)');
ylabel('Height (cm');
title('Direct Numerical Solution: Temperature (C)');
caxis([0 100]);
colorbar;

%Apply 2D PDE solution, iteratively, via Jacobi Iteration
n=length(x);
m=length(y);
Ti=zeros(size(X));
Ti(end,:)=T0*sin(pi*x/(w));
T=zeros(size(X));
toler=1e-2;
iter=0;
while abs(norm(Ti)-norm(T))>toler
    T=Ti;
    iter=iter+1;
    for i=2:m-1
        for j=2:n-1
            Ti(i,j)=T(i,j)+(1/4)*(T(i-1,j)+T(i+1,j)+T(i,j-1)+T(i,j+1)-4*T(i,j));
        end;
    end;
end;

Ti
iter

%Plot Results
figure(3);
imagesc(x,y,Ti); 
axis xy;
xlabel('Width (cm)');
ylabel('Height (cm');
title('Iterative Numerical Solution: Temperature (C)');
colorbar;
