function transient_sq_conduction_4_T3 (bc_value, ic_value)
%  1/8 symmetry of sq with heat generation Q          6
% connections: (1)  1  2  3                         /  |
%              (2)  2  4  5                       3----5
%              (3)  3  5  6                     / |  / |
%              (4)  5  3  2                   1---2/---4

%          part data:
n_side = 2          ; % number of elements per side
L      = 4 / n_side ; % element length in x & y
rho_cp = 1          ; % material density * specific heat
k      = 8          ; % material isotropic thermal conductivity
Q      = 6          ; % heat generation rate per unit volume
thick  = 1          ; % part thickness
n_bc   = 3          ; % number of essentrial BC
n_unkn = 3          ; % number of unknowns
R_bc   = ones(n_bc,1) * bc_value ; % known BC temperatures

%   generalized trapezoidal time integration data
beta   = 0.5 ; % Crank-Nicolson, Galerkin=2/3,
               % foward difference=0, backward difference=1
n_steps = 25  ; % number of times steps
delta_t = 0.1 ; % size of time step
time = [0:1:n_steps] * delta_t  ; % actual time values
R_u = zeros (n_unkn, n_steps+1) ; % allocate time history
R_u (:, 1) = ic_value           ; % insert initial condition

%  partitioned assembly, with unknown values first
%  M_uu * R_u_dot + S_uu * R_u = C_Q + C_bc
S_uu = [ 1,  -1,  0 ; ...
        -1,   4, -2 ; ...
         0,  -2,  4 ] * k * thick / 2           ; % conduction
M_uu = [ 2,   1,  1 ; ...
         1,   6,  2 ; ...
         1,   2,  6 ] * rho_cp * thick * L^2 / 24 ; % capacity
C_Q  = [ 1,   3,  3 ]' * Q * thick *L^2 / 6 ; % internal source

%  known partition for BC contribution
S_ub = [ 0,  0, 0 ; ...
        -1,  0, 0 ; ...
         0, -2, 0 ] * k * thick / 2 ; % conduction
C_bc = -S_ub * R_bc  ; % move BC conribution to RHS

%  form matrices for linear system S_hat * R_u(time) = F(time)
S_hat = M_uu / delta_t + beta * S_uu  ; % constant in time
S_inv = inv (S_hat)    ; % invert only once, lu is better here
S_RHS = S_hat - S_uu   ; % constant in time

%  time marching solution, allow Q to vary in time
for j = 2:1:n_steps+1         ; % march to next time
  scale = 1 + 0 * time (j)^2  ; % fake time history of Q
  F = C_Q * scale + C_bc + S_RHS * R_u (:, j-1) ; % at time
  R_u (:, j) = S_inv * F                        ; % answer
% React (:, j) = S_ub' * R_u (:, j) + S_bb * R_bc ...
%              - C_Q * scale ; % optional reaction time history
end % for all time steps

%  plot selected time histories
clf  % clear figure
hold on
title  ('First Three Nodal Time Histories')
xlabel ('Time (sec)')
ylabel ('Temperature (F)')
plot ( time, R_u(1,:),'ko') ; plot ( time, R_u(2,:),'kd')
plot ( time, R_u(3,:),'ks')
legend ('Node 1', 'Node 2', 'Node 3')
plot ( time, R_u(1,:),'k-') ; plot ( time, R_u(2,:),'k-')
plot ( time, R_u(3,:),'k-')
grid on
% end transient_sq_conduction_4_T3