% Example 13.6 (Unsteady) Heat Equation
% Semi-discretization Method

clear all
close all
clc

% Physical parameters
alpha = 1;
L = 1;

% Mesh
N=21;
dx = 1/(N-1);
x = 0:dx:L;

% time advancement
dt = 0.00125;
tmax = 2;
t = 0:dt:tmax;
nmax = length(t);

T = zeros(nmax,N); % Allocate solution matrix

% initial condition
T(1,:) = sin(pi().*x);

% inhomogeneous part
fvec = @(x,t) (pi().^2-1).*exp(-t).*sin(pi().*x);

% Matrix
twos = -2*ones(1,N-2);
A = diag(twos)+diag(ones(1,N-3),-1)+diag(ones(1,N-3),1);

%Explicit Euler
for n = 1:nmax
    dTdt = (alpha/dx^2).*A*T(n,2:N-1)' + fvec(x(2:N-1),t(n))';
    T(n+1,2:N-1)=T(n,2:N-1)' + dt.*dTdt;    
end

figure
%hold on
for n = 1:nmax
    plot(x,T(n,:),'-k')
    ylim([0,1])
    xlabel('x')
    ylabel('T')
    title(['Explicit Euler solution of Heat Equation with dt = ',num2str(dt)])
    %saveas(gcf,[num2str(n),'.png'])
    pause(dt)
end

tplots = [1, 401, 801, 1201, 1601];
figure
hold on
for n = 1:length(tplots)
    plot(x,T(tplots(n),:),'LineWidth',2)     
end
ylim([0,1])
xlabel('x')
ylabel('u')
title(['Explicit Euler solution of Heat Equation with dt = ',num2str(dt)])
legend(['t = ',num2str(t(1))],['t = ',num2str(t(401))],['t = ',num2str(t(801))],['t = ',num2str(t(1201))],['t = ',num2str(t(1601))])






