function Uxx_U_Q_L3_GenBC (n_e, show)
% File: Uxx_U_Q_L3_GenBC.m for E*Uxx + k*U + Q(x) = 0, E=1=k, plus
% two general EBC or 1 EBC and 1 NBC.  SEE COMMENTS AT THE BOTTOM.
% if U(0)=0=U(L) then U(x) = sin(x)/sin(L) - x
% if U(0)=0, dU/dx(L)=0 then U(x) = sin(x)/cos(L) - x
% if U(0)=0, dU/dx(L)=-0.357907 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   x_5= L
%  U_o  *------(1)-----*------(2)------*  U(L)
%       T_1    T_2    T_3     T_4     T_5
% Connectivity list:    e     i      j     k
%                       1     1      2     3  % L3
%                       2     3      4     5  % L3
%.....................................................................
if ( nargin == 1 ) ; 
  show = 0 ;             % flag to show exact solution (if = 1)
elseif ( nargin == 0 ) ; 
  show = 0 ; n_e = 2 ;   % n_e = number of elements
end % if from argument count

%            Set given constants (try changing n_e)
n_g = 1 ;                  % number of DOF per node
n_n = 3 ;                  % number of nodes per element
n_i = n_n*n_g           ;  % number of DOF per element
n_m = 1 + n_e*(n_n - 1) ;  % number of mesh nodes
n_d = n_g*n_m           ;  % system degrees of freedom (DOF)

% specific problem data
L = 1.0 ; L_e = L/n_e               ; % domain and element lengths
E = 1.0 ; k = 1.0                   ; % stiffness data
X = [0:1:(n_d-1)]*L/(n_d-1)         ; % sequential element lengths

%   Allocate (wasteful) storage for BC flags and values
EBC_flag  = zeros(n_m,n_g) ;          NBC_flag  = zeros(n_m,n_g) ; 
EBC_value = zeros(n_m,n_g) ;          NBC_value = zeros(n_m,n_g) ;

%   Manually set specific essential & Neumann bc's (SEE BOTTOM)
EBC_flag (1,   1) = 1 ;            EBC_value (1,   1) = 0.       ;
EBC_flag (n_m, 1) = 1 ;            EBC_value (n_m, 1) = 0.       ;
% NBC_flag (n_m, 1) = 1 ;            NBC_value (n_m, 1) = 0.       ;

%         Constant quadratic element square matrices
K_e = E   * [ 7, -8,  1  ; ...              % axial
             -8, 16, -8  ; ...
              1, -8,  7 ] / (3*L_e)  ;
M_e = L_e * [4, 2, -1 ; ...                 % mass
             2, 16, 2 ; ...
            -1,  2, 4 ] / 30         ;
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
  rows (1:n_i) = [1:n_i] + (n_i-1)*(j-1)   ; % element DOF numbers
  S (rows, rows) = S (rows, rows) + S_e    ; % add to system sq
  Q_e = [X(rows)]'                         ; % nodal source Q(x)=x
  C_e = M_e * Q_e                          ; % element source
  C (rows) = C (rows) + C_e                ; % add to sys column
end % for

% Optional storage allocation for reaction recovery
EBC_count = sum(sum(EBC_flag))                  ; % number of EBC's
NBC_count = sum(sum(NBC_flag))                  ; % number of NBC's
EBC_row   = zeros(EBC_count, n_d) ; EBC_col = zeros(EBC_count, 1) ;
EBC_react = zeros(EBC_count, 1)   ;     % reaction arrays allocated

%  Apply essential BC at nodes, via trick to avoid partitions
EBC_flag  = reshape (EBC_flag,  n_d, 1)        ; % change to vector 
NBC_flag  = reshape (NBC_flag,  n_d, 1)        ; % change to vector 
EBC_value = reshape (EBC_value, n_d, 1)        ; % change to vector 
NBC_value = reshape (NBC_value, n_d, 1)        ; % change to vector 
kount = 0                                      ; % initialize counter
for j = 1:n_d       % check all DOF for essential and natural BC
  if ( NBC_flag(j) )                             % then NBC here
    C (j) = C (j) + NBC_value (j)              ; % add NBC value
  end % if NBC for this DOF

  if ( EBC_flag(j) )                         % then EBC here
    if ( EBC_flag(j) == NBC_flag(j) )        % warning both here
      fprintf ('Usually fatal error: EBC and NBC at same DOF \n')
    end % if fatal BC data

%   Save reaction data to be destroyed by EBC solver trick
    kount = kount + 1                      ; % copy counter
    EBC_row(kount, 1:n_d) = S (j, 1:n_d)   ; % copy reaction data
    EBC_col(kount)        = C (j)          ; % copy reaction data

%   Insert EBC identity to avoid matrix partition accounting
    EBC = EBC_value(j)                ; % recover EBC value
    C (:) = C (:) - EBC * S (:, j)    ; % carry known column to RHS
    S (:, j) = 0 ; S (j, :) = 0       ; % clear, restore symmetry
    S (j, j) = 1 ; C (j) = EBC        ; % insert identity into row
  end % if EBC for this DOF
end % for over all DOF

T = S \ C                           % Solve for unknown values
EBC_react = EBC_row * T - EBC_col   % Reaction from matrix system (+ or -)

% End finite element calculations. ---- Begin graphics ----
 
clf ; hold on ; grid on                         % clear & keep
title ('E*Uxx + k*U + Q(x), E=k=1, Q(x) = x')

%  Plot the exact results
z = [0:0.025:1] * L ;                           % points for graph
if (EBC_count == 2 )                            % for 2 EBC
  exact = sin(z)/sin(L) - z                   ; % analytic
else  % 1 EBC, 1 NBC                            % 1 EBC 1 NBC
  exact = sin(z)/cos(L) - z                   ; % analytic
end % if true solution
plot (z, exact, 'r.-');  plot (X, T, 'ko')      % exact & nodal
legend ('Exact', 'FEA')                       ; % legend

%        Show BC locations on the plot
k_EBC = 0 ; k_NBC = 0 ;
for j = 1:n_d % check for essential and natural BC
  node  = 1 + (j - 1)/n_g                     ; % node at DOF j
  first = n_g*(node - 1) + 1                  ; % first DOF at node
  if ( EBC_flag(j) )                            % then EBC here
    k_EBC = k_EBC + 1                         ; % counter inc
    b_text = sprintf('--EBC (%g), React: %g', k_EBC, -EBC_react(k_EBC));
    plot (X(node), T(first), 'k*')              % star at EBC node
    text (X(node), T(first), b_text)            % notes at EBC node
  end % if                                      % EBC completed   
  if ( NBC_flag(j) )                            % then NBC here
    k_NBC = k_NBC + 1                         ; % counter inc
    b_text = sprintf('--NBC (%g) = g', k_NBC, NBC_value(j)) ; 
    plot (X(node), T(first), 'k*')              % star at NBC node
    text (X(node), T(first), b_text)            % notes at NBC node
  end % if NBC at this node
end % for over all DOF

%  ------------------------------------------------------
%  graph of value in 3 node quadratic line elements
%  ------------------------------------------------------
% e_x, e_y   = coordinates & data values at 3 element nodes
% x_el, y_el = interpolated parametric values

 n_poly = ceil (75/n_e) ; % curve segments
 e_x(3)=0 ; e_y(3)=0 ; H(3)=0 ; % pre-allocation

%    Loop over all L3 elements
  V_X = -realmax ; V_N = realmax ; % for interpolated values
  for it = 1:n_e  ;
    
%      Extract nodal coordinates & values. Plot nodes.
    rows = [1, 2, 3] + (it - 1)*2     ;  % get connectivity 
    e_x  = X (rows) ; e_y  = T (rows) ;  % gather element data
    plot (e_x, e_y, 'ko')                % nodal symbols
    e_text = sprintf ('   (%g)', it)   ; % offset # from pt
    text (e_x(2), e_y(2), e_text)        % Plot element number

%     Loop over local points on the quadratic element
    for k = 1: (n_poly + 1) % points in parametric space
%     get parametric interpolation functions
      Z    = 2*(k - 1)/n_poly - 1            ; % on -1 to 1
      H(1) = 0.5*(Z*Z - Z) ; H(2) = 1. - Z*Z ;
      H(3) = 0.5*(Z*Z + Z) ;
      x_el (k) = dot(H,e_x) ; y_el (k) = dot(H,e_y) ; % true value
    end % for k parametric points in element
    plot (x_el, y_el, 'k-') % parametric curve through 3 nodes
    if ( max(y_el) > V_X )
      V_X = max(y_el) ;     % curve exceeds max node value
    end % if
    if ( min(y_el) < V_N )
      V_N = min(y_el) ;     % curve below min node value
    end % if
  end % for it over all elements
  xlabel (['Position (of ', int2str(n_e), ' elements)']) ; 
  ylabel (['U(x). FE Min, Max = ', ...
            num2str(V_N), ', ', num2str(V_X)])
% end of Uxx_U_Q_L3_GenBC

% ------------- Begin Samples --------------------------------
% if U(0)=0=U(L) then U(x) = sin(x)/sin(L) - x
%  EBC_flag (1,   1) = 1 ;      EBC_value (1,   1) = 0. ;
%  EBC_flag (n_m, 1) = 1 ;      EBC_value (n_m, 1) = 0. ;

% n_e = 1, EBC_react = -0.1898 -0.3565
% T = 0 0.0694 0

% n_e = 3, EBC_react = -0.1884 -0.3579
% T = 0 0.0305 0.0555 0.0697 0.0682 0.0463 0

% if U(0)=0, dU/dx(L)=0 then U(x) = sin(x)/cos(L) - x
%  EBC_flag (1,   1) = 1 ;      EBC_value (1,   1) = 0. ;
%  NBC_flag (n_m, 1) = 1 ;      NBC_value (n_m, 1) = 0. ;

% n_e = 1, EBC_react = -0.8489 
% T = 0 0.3849 0.5540
   
% n_e = 2, EBC_react = -0.1882
% T = 0 0.0439 0.0696 0.0598 -0.0002

% if U(0)=0, dU/dx(L)=-0.357907 then U(x) = sin(x)/sin(L) - x
%  EBC_flag (1,   1) = 1 ;      EBC_value (1,   1) = 0.       ;
%  NBC_flag (n_m, 1) = 1 ;      NBC_value (n_m, 1) = -0.35797 ;

%  n_e = 2, EBC_react = -0.1882
%  T = 0 0.0439 0.0696 0.0598 -0.0002

%  n_e = 4, EBC_react = -0.1883
%  T = 0 0.0231 0.0440 0.0602 0.0697 0.0703 0.0600 0.0370 -0.0001
% -------------- end samples -----------------------------------