%FInite difference method to solve the heat transfer problem in chapter 8
clc;    clear all;

T_0 = [80;30;30;30;30];    %Initial temperature

k = 230; %input('Thermal conductivity of the rod (in W/m C): '); 
C = 900; %input('Specific heat of the rod (in J/Kg C): '); 
rho = 2700; %input('Density of the rod (in Kg/m^3): '); 

t = 100;    %input('Total time for which temperature profile is required (in s): ');
del_t = 0.1; %input('Time interval or increment of march in time axis (in s): ');
del_x = 0.025; %input('Length of each element (in m): ');

iterations = t/del_t;  


T = zeros(iterations+1,5);   %Temperature matrix, column is different time step, row is the temperature value at each time step

D = k/(C*rho);     %diffusivity constant

%-------Initial temperaure condition------------------%
for i = 1:5
    T(1,i)=T_0(i);
end

%-----solution or time march---------------%
for i =1:iterations
    T(i+1,1) = 80;         %fixed temperatures
    T(i+1,5)= 30;
    
    for j = 2:4
        T(i+1,j)=T(i,j)+(D*del_t/del_x^2)*(T(i,j-1)-2*T(i,j)+T(i,j+1));    %value of temperature using finite different analysis
    end
end

plot(T);