function [] = Example_Transients () % revised 4/20/20
% Time stepping example of a transient conducting bar
%   with a mesh of three L2 linear line elements
% PDE:   K A d2T/dx2 - Q = rho cp A dT/dt
% Matrix: [M]{T_dot} + [S]{T} = c_Q(t)
%
% Non-dimensional form [S*] + a [M*] = {c*(t)}
%             a = rho cp Le^2 / dt K A
%     See text Fig 15.3 for time history graphs
%----------------------------------------------------
% EBC:  V_0 1         2         3         4 V_L
%           .***(1)***.***(2)***.***(3)***.      Q=0
% ICs:     T_0       T_0       T_0       T_0
%              Only T_2 and T_3 are free
%------------------------MESH------------------------
I_C = 0 ; T_0 = I_C*[1; 1]     ; % Initial Conditions
V_0 = 10 ; V_L = 20           ; % Constant EBC values
x=[0 1 2 3]'                        ; % x Coordinates
n_steps = 30                 ; % Number of time steps
T = zeros (2, n_steps) ; % Allocate Free temperatures
diag = 1 ;% turn on/off use only diagonal mass matrix
ave  = 0 ;% turn on/off use only averaged mass matrix
a = 8.0              ; % *** Bigger for smaller Del_T

K = [2   -1; -1  2]   ; % Assembled conduction matrix
M = [4  1;  1  4]*a/6   ; % Assembled capacity matrix
Q = [0 ; 0]               ; % Assembled source vector

if ( diag == 1 )           ; % then use diagonal mass
   M = [5  0;  0  5]*a/6   ; % scaled diagonal of {M}
end                       ; % if diagonal mass matrix

if ( ave == 1 )            ; % then use averaged mass
  M = ([4  1;  1  4]*a/6 + [5  0;  0  5]*a/6)/2;% ave
end                       ; % if averaged mass matrix

S = K + M          ; % New linear system at each time
S_inv = inv(S)  ; % Invert only once (constant Del_T)
BC_1 = V_0*[-1;0]; BC_4 = V_L*[0;-1]; % c_EBC vectors
c = M*T_0 - BC_1 - BC_4          ; % first c(t=Del_T)

T(:,1) = S_inv * c ; % solve for first history result

for n = 2:n_steps ; % time step loops ==> ==> ==> ==>
  c = M*T(:, n-1)-BC_1-BC_4 ; % update the RHS source
  T(:,n) = S_inv * c ; % solve for T at this new time
end ; % for each step after the first <== <== <== <==
% ----------------end time marching -----------------

%  Restore EBC to graph mesh Show values through time
Show = zeros (n_steps+1, 4)   ; % allocate graph data
Show(1,:) = I_C   ; % place ICs in the first row, t=0
Show(:,1) = V_0    ; % put left BC in col 1 for all t 
Show(:,4) = V_L;  ; % put right BC in col 4 for all t
Show(2:n_steps+1, 2:3) = T'  ; % insert computed T(t)
%  Scale the axes, label axes, set title
xmin = 0; xmax = 3;          % set length axis bounds
ymin = min([T_0; V_0; V_L])   ; % set min temperature
ymax = max([T_0; V_0; V_L])   ; % set max temperature
clf; hold on ; grid on ; axis([xmin xmax ymin ymax]);
xlabel('x'); ylabel('Temperature')     ; % label axes
title('FEA Transient Conduction Solution')  ; % title
%   Add stready state temperatures for reference
end_x=[0 3];end_t=[V_0 V_L];% steady state end values
plot(end_x, end_t,'k--')    ; % show dashed reference
%   Graph each t state approach steady state value
for n = 1:n_steps+1; % loop over time --> --> --> -->
  plot (x, Show(n,:),'k-')    ; % show element values
  plot (x, Show(n,:),'k*')       ; % show node values
end; % for inital and following steps <-- <-- <-- <--
plot (x, Show(1,:),'ko'); % repeat initial conditions
ave=(V_0+V_L)/2                     ; % text location
text(0.5, ave,['  IC = ', num2str(I_C)])  ; % list IC
print -dpng bar_transients; % save to a png plot file
fprintf ('Created file bar_transients.png \n')
hold off                                 ; % all done
% end Example_Transients ============================