%% 1D Parabolic PDE - Steady-State, Explicit, and Implicit Heat Eqn.
% Solving dT/dt=alpha*Txx

%Domain in Centimeters
w=15;
T0=100;
dx=1.0;
x=[0:dx:w];
Ta=x.*(T0/w)   %Linear increase from one side to other

%% Steady-State (Elliptic) solution to BVP, T=0 at x=0, T=100 at x=w
n=length(x)-2;
C=zeros(1,n)'; C(end)=-T0;
diagonal=sparse(1:n,1:n,-2*ones(1,n),n,n);
suprdiag=sparse(2:n,1:n-1,ones(1,n-1),n,n);
smatrix=suprdiag+diagonal+suprdiag';
fmatrix=full(smatrix);
Tss=fmatrix\C;
Tss=[0; Tss; T0]';  %Set boundaries
Tss

figure(1);
plot(x,Tss);
xlabel('Width (cm)');
ylabel('Temperature (C)');
title('Steady-State Solution');

%% Explicit (Parabolic) solution to Heat Eqn., T=0 at x=0, T=100 at x=w
m=length(x);
kappa=1.11;  %Conductivity of copper in cm^2/s
dt=1.0*dx^2/(2*kappa);
t=0;
toler=1;
T1=zeros(size(x));
T1(end)=T0;
T=zeros(size(x));
while abs(norm(T1)-norm(Ta))>toler
    T=T1;
    t=t+dt;
  %  for i=2:m-1
  %      T1(i)=T(i)+(kappa*dt/(dx^2))*(T(i-1)-2*T(i)+T(i+1));
  %  end;
    T1(2:m-1)=T(2:m-1)+(kappa*dt/(dx^2)).*(T(1:m-1-1)-2.*T(2:m-1)+T(3:m-1+1));
    figure(2);
    plot(x,Ta,x,T1);
    legend('Analytical S-S Solution','Numerical Solution');
    xlabel('Width (cm)');
    ylabel('Temperature (C)');
    title(['Explicit Solution, t=',num2str(t)]);
    pause(0.001);
end;
T1
t

%% Implicit (Parabolic) solution to Heat Eqn., T=0 at x=0, T=100 at x=w
m=length(x)-2;
kappa=1.11;  %Conductivity of copper in cm^2/s
dt=10*dx^2/(2*kappa);
t=0;
toler=1;
T1=zeros(size(x))';
T1(end)=T0;
T=zeros(size(x))';
while abs(norm(T1)-norm(Ta))>toler
    T=T1;
    t=t+dt;
    d=kappa*dt/(dx^2);
    diagonal=sparse(1:m,1:m,(1+2*d)*ones(1,m),m,m);
    suprdiag=sparse(2:m,1:m-1,-d*ones(1,m-1),m,m);
    smatrix=suprdiag+diagonal+suprdiag';
    fmatrix=full(smatrix);
    T(2)=T(2)+T(1)*d;
    T(end-1)=T(end-1)+T(end)*d;
    T2=fmatrix\T(2:end-1);
    T1=[0; T2; T0];
 %  T1(1)=0; T1(end)=T0;
    figure(3);
    plot(x,Ta,x,T1);
    legend('Analytical S-S Solution','Numerical Solution');
    xlabel('Width (cm)');
    ylabel('Temperature (C)');
    title(['Implicit Solution, t=',num2str(t)]);
    pause(0.001);
end;
T1
t
        
 
%% Implicit Crank-Nicolson solution to Heat Eqn., T=0 at x=0, T=100 at x=w0
m=length(x)-2;
kappa=1.11;  %Conductivity of copper in cm^2/s
dt=10*dx^2/(2*kappa);
t=0;
toler=1;
T1=zeros(size(x));
T1(end)=T0;
T=zeros(size(x));
while abs(norm(T1)-norm(Ta))>toler
    T=T1;
    t=t+dt;
    d=kappa*dt/(dx^2);
    for i=2:(m+1)
       T1(i)=d*T(i-1)+2*(1-d)*T(i)+d*T(i+1);
    end;
    T1(2)=T1(2)+T(1)*d;
    T1(end-1)=T1(end-1)+T(end)*d;
    diagonal=sparse(1:m,1:m,2*(1+d)*ones(1,m),m,m);
    suprdiag=sparse(2:m,1:m-1,-(d)*ones(1,m-1),m,m);
    smatrix=suprdiag+diagonal+suprdiag';
    fmatrix=full(smatrix);
    T2=(fmatrix\T1(2:end-1)');
    T1=[0; T2; T0]';  %Set boundaries
    figure(4);
    plot(x,Ta,x,T1);
    legend('Analytical S-S Solution','Numerical Solution');
    xlabel('Width (cm)');
    ylabel('Temperature (C)');
    title(['Implicit Solution, t=',num2str(t)]);
    pause(0.001);
end;
T1
t


        