%Finite element method to solve the transient heat transfer problem in chapter 8
%only 4 elements are allowed; mainly to see the effect of time march and
%the initial conditions on how fast the solution converges

clc;    clear all;

%---------inputs----------%
T1 = input('Temperature at the left end of the rod (in deg C): ');
T_before = input('Temperature of the rod before heating started (in deg C): ');

t = input('Total time for which temperature profile is required (in s): ');
del_t = input('Time interval or increment of march in time axis (in s): ');

%----------calculations------------------------%
iterations = t/del_t;      %total number of iterations
T_0 = [T1;T_before;T_before;T_before;T_before];   %Initial temperature

T = zeros(iterations+1,5);    %temperature matrix, columns are time steps, rows are the temperature of 5 nodes at a time step
time = zeros(iterations+1,1);

%dummy storage vectors for storing active components are temperatures
T_d1 = zeros(3,1);            
T_d = zeros(3,1);


for i = 1:5                    %Initiatal temperature
    T(1,i)=T_0(i);
end

A=[0.498 -0.3521 0.0880; -0.3521 0.5869 -0.3521; 0.0880 -0.3521 0.4988];   %see lecture notes for its definition; 
f = [18.0456;-6.4558;7.7762];                                              %see lecture notes for its definition; 

for i =2:iterations+1
    for j = 1:3
        T_d1(j) = T(i-1,j+1);    %value of temperature from last time step
    end
    
    T_d = T_d1 - A*T_d1*del_t+f*del_t;   %calculation of temperature at the current time step
    
    T(i,1) = 80;           %fixed temperature at left end
    T(i,5) = 30;           %fixed temperature at right end
    
    for j = 1:3            %temperature of "active" nodes at the current time step
        T(i,j+1)=T_d(j);
    end
    time(i)= i*del_t;      %time step
end

%Plotting temperature profile with time
plot(time, T);
axis([0 t T_before-10 T1+10]);
xlabel('Time');
ylabel('Temperature (deg C)');