
% Program to illustrate various FD approximations 
% for the Cooling problem

% remove all variables and close all windows
clear
close all


T0=1.                           % Set initial temperature
tau=.7;                         % Set thermal conductivity
dt=input(' Give dt : ');        % Get time step interactively
tmax=10;                        % Give maximum time to calculate

nt=round(tmax/dt);              % Calculate required number of time steps 

% initialize vectors needed for calculations
% T -> simplest FD
% Ta -> analytical solution
% Ti -> implicit
% Tm -> mixed
T=zeros(nt,1);
Ta=T;
Ti=T;
Tm=T;

% set initial temperature for all methods
T(1)=T0;
Ta(1)=1;
Ti(1)=1;
Tm(1)=1;

for i=1:nt,
t(i)=i*dt;
T(i+1)=T(i)-dt/tau*T(i);                     % explicit forward
Ta(i+1)=exp(-dt*i/tau);                      % analytic solution
Ti(i+1)=T(i)*(1+dt/tau)^(-1);                % implicit 
Tm(i+1)=(1-dt/(2*tau))/(1+dt/(2*tau))*Tm(i); % mixed
end

plot(t,T(1:nt),'b-',t,Ta(1:nt),'r-',t,Ti(1:nt),'g-',t,Tm(1:nt),'k-')
xlabel('Time(s)')
ylabel('Temperature')
title(sprintf('dt=%g; tau=%g',dt,tau))