%% Separation of variables solution of heated bar problem
%% Left end at Tin; right end at 0; initial temperature = To
%% Demonstrating linear superposition of solutions
Nmodes = 40;
Nx = 41;
Nt = 21;
tmax = 0.2;

Tin = 100;
To = 20;

n = 1:Nmodes;
x = linspace(0, 1, Nx);
t = linspace(0, tmax, Nt);
Xmode = sin(n' * pi * x);   
Tpart = exp(-(n' * pi).^2 * t)';

figure(1), clf

%% Steady state solution
Uss = ones(Nt, 1) * Tin * (1 - x);
subplot(1,3,1)
surf(x, t, Uss);
xlabel('x / L'), ylabel('(D/L^2) t'), zlabel('Temp'), title('STEADY-STATE')
axis square

%% Solution with homogeneous BCs and modified ICs
bn = To* 2/pi * diag((1 - (-1).^n)./n) - Tin * 2/pi * diag(1./n) ;
Ubc = Tpart * bn * Xmode;
subplot(1,3,2)
surf(x, t, Ubc)
xlabel('x / L'), ylabel('(D/L^2) t'), zlabel('Temp'), title('TRANSIENT')
axis square

%% Total solution = steady-state + transient
U = Uss + Ubc;
subplot(1,3,3)
surf(x, t, U)
xlabel('x / L'), ylabel('(D/L^2) t'), zlabel('Temp'), title('TOTAL SOLUTION')
axis square

figure(2), clf
%% Total solution plotted as time slices
tslices = [1, 2, 5, 10, 21];
plot(x, U(tslices, :))
xlabel('x / L')
ylabel('Temperature')
legtext = [char(ones(size(tslices))'*'(D/L^2) t = '), strvcat( num2str( t(tslices)' ) ) ];
legend(legtext)
title('HEATED BAR')