% File: Uxx_U_Q_2EBC.m for ODE E*Uxx + k*U + Q(x) = 0, E=1=k, plus 2 EBC
% U(0)=0=U(L) then U(x) = sin(x)/sin(L) - x
% Ref: J.E. Akin, "FEA w Error Estimators" Elsevier, 2005, page 50  
%.......................................................................
%      x_1=0     x_2       x_3       x_4 = L
%  U_o  *---(1)---*---(2)---*---(3)---*  U(L)
%       T_1       T_2       T_3       T_4
% Connectivity list:    e     i      j
%                       1     1      2   % L2
%                       2     2      3   % L2
%                       3     3      4   % L2
%.......................................................................
%            Set given constants (try changing n_e)
n_e = 3 ; n_d = n_e + 1 ;  % number of elements and nodes (DOF)
n_n = 2                 ;  % number of nodes per element
L = 1.0 ; L_e = L/n_e   ;  % domain and element lengths
E = 1.0 ; k = 1.0       ;  % stiffness data
U_o  = 0.0 ; U_L  = 0.0 ;  % essential boundary conditions
X = [0:n_e] * L_e       ;  % merge lengths

%         Constant element square matrices
K_E = E   * reshape ([1, -1, -1, 1], 2, 2) / L_e; % axial
M_e = L_e * reshape ([2,  1,  1, 2], 2, 2) / 6  ; % mass
K_k = k * M_e                                   ; % foundation
S_e = K_E - K_k                                 ; % net square matrix
Q_e = zeros (n_n,1) ; C_e = zeros (n_n,1)       ; % pt & elem sources

%         Assemble n_d by n_d square matrix terms
S = zeros (n_d, n_d) ;  C = zeros (n_d, 1) ; % initalize sums
for j = 1:n_e                          ; % loop over elements
  last = j + n_n - 1 ; first = last - n_n + 1 ; % DOF numbers
  rows = [first last]                  ; % vector subscripts
  S (rows, rows) = S (rows, rows) + S_e; % add to system sq
  Q_e = [X(first);  X(last)]           ; % nodal source
  C_e = M_e * Q_e                      ; % element source
  C (rows) = C (rows) + C_e            ; % add to sys column
end % for

%  Save reaction data to be destroyed by EBC solver trick
Row_1_S   = S(1,   1:n_d) ; Row_1_C   = C (1) ;
Row_n_d_S = S(n_d, 1:n_d) ; Row_n_d_C = C (n_d) ;

%  Apply essential BC at node 1, via trick to avoid partitions
C (:) = C (:) - U_o * S (:, 1)   ; % carry known column to RHS
S (:, 1) = 0 ; S (1, :) = 0      ; % clear, restore symmetry
S (1, 1) = 1 ; C (1) = U_o       ; % trick not to change array size
C (:) = C (:) - U_L * S (:, n_d) ; % carry known column to RHS
S (:, n_d) = 0 ; S (n_d, :) = 0  ; % clear, restore symmetry
S (n_d, n_d) = 1 ; C (n_d) = U_L ; % trick not to change array size
T = S \ C                        ; % Solve for four values
q_o = -( Row_1_S * T - Row_1_C ) ; % Get reaction from matrix system
q_L =  ( Row_n_d_S * T - Row_n_d_C ); % reaction from matrix system
 
%            Plot the temperature results
clf ; hold on ; grid on                    % clear & keep
title ('ODE E*Uxx + k*U + Q(x), E=k=1, Q(x) = x')
text  (mean(X), mean(T), ' U(0)=0,   U(L)=0 ')
q_text = ['Reactions: ', num2str(q_o), ', and ', num2str(q_L)];
text  (mean(X)/2, mean(T)/2, q_text) 
xlabel ('Position') ; ylabel ('U(x)')
x = [0:0.025:1] * L ; exact = sin(x)/sin(L) - x ; % analytic
plot (x, exact, 'b.-');  plot (X, T, 'ko')       % exact & nodal
legend ('Exact', 'FEA'); plot (X, T, 'k-')       % legend & elements
% end