% Single quadratic line element for conducting-convecting rod
%              (e) k,A,h,p,T_inf
% T_0  *---------*--------* insulated
%      1 x=0     2        3 x=L
L = 4; A = 0.01389; h = 2 ; k = 120;      % data
p = 0.5; T_0 = 10; T_inf = 0;             % data
x_e = [0, 2, 4];                          % node coordinates
a = k*A/3/L ; b = h*p*L/30 ;              % shorthand
%         System Matrices
K_e = [7, -8, 1; -8, 16, -8; 1, -8, 7]*a  % conduction
%	K_e = 0.9723   -1.1112    0.1389
% 	     -1.1112    2.2224   -1.1112
% 	      0.1389   -1.1112    0.9723
M_e = [4, 2, -1; 2, 16, 2; -1, 2, 4]*b    % convection
%	M_e = 0.5333    0.2667   -0.1333
%          0.2667    2.1333    0.2667
%         -0.1333    0.2667    0.5333
C_e = 5*b*T_inf*[1; 4; 1]                 % source
%	C_e' = 0 0 0
K = K_e + M_e                             % assembly
%    K = 1.5056   -0.8445    0.0056
%       -0.8445    4.3557   -0.8445
%        0.0056   -0.8445    1.5056
%         Apply BC
C = C_e - T_0*K(:,1)                      % BC to RHS
%    C' = -15.0563 8.4453 -0.0557
%         Solve
D(2:3) = K(2:3,2:3)\C(2:3)          % solve partition
D(1) = T_0                          % BC identity
%    D' = 10.0000    2.1675    1.1788
%         Reaction heat flux (with C before BC)
q = K*D' - C_e
%    q' = 13.2324 0.0000 0.0000     % heat flow in
%         Post-process
Loss = 5*b*[1, 4, 1]*D'             % total, exact = 12.9
%    Loss = 13.2324                 % heat flow out
%         Compare
x =  [0:0.10:4];                    % domain
a = [-1:0.05:1];                    % parametric
approx = 0.5*(a.^2 - a)*D(1) + (1 - a.^2)*D(2) ...
       + 0.5*(a.^2 + a)*D(3) ;      % interpolate
m = sqrt(h*p/(k*A));                % for exact
exact = T_inf + (T_0-T_inf)*cosh(m*(L-x))/cosh(m*L);
%         Graph
clf                                 % clear graphics
plot(x, approx,'k-')                % approximate
hold on
plot(x,exact,'b--')                 % graph exact
plot(x_e, D,'ko')                   % graph node values
legend ('Single element', 'Exact solution')
grid                                % reference grid
xlabel ('Position (ft)') ; ylabel ('Temperature (F)')
title ('Axially Conducting and Convecting Bar')