% The function illustrates midpoint-based second order 
% accurate finite elements;

function fem_02
  
% Size of the problem;
N1 = 6;

% Viewing the basis functions
Nx = 200;
dx = 1.0 / (Nx - 1);
X = 0.0 : dx : 1.0;

figure(1); 
plot( X, phi( N1, 1, X ), X, phi( N1, 2, X ), X, phi( N1, 3, X ), ...
      X, phi( N1, N1, X ) ); hold on;
plot( X, psi( N1, 1, X ),':', X, psi( N1, 2, X ),':', X, psi( N1, 3, X ),':', ...
      X, psi( N1, N1, X ),':' ); hold off;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



% The basis functions based on the interval endpoints;
function [fun] = phi( N, k, x )
  
  Nx = length( x );
  
  h = 1.0 / N;
  for i=1:Nx
    if ( (x(i) > (k-1) * h) & (x(i) <= k * h) )
      xi = ( k * h + (k-1) * h ) / 2.0;
      fun(i) = 2.0 * ( x(i) - (k-1) * h ) * ( x(i) - xi ) / h^2;
    elseif ( (x(i) > k * h) & (x(i) <= (k+1) * h) )
      xi = ( (k+1) * h + k * h ) / 2.0;
      fun(i) = 2.0 * ( x(i) - (k+1) * h ) * ( x(i) - xi ) / h^2;
    else
      fun(i) = 0.0;
    end;
  end;
  
  
% The basis functions based on the midpoints;
function [fun] = psi( N, k, x )
  
  Nx = length( x );
  
  h = 1.0 / N;
  for i=1:Nx
    if ( (x(i) > (k-1) * h) & (x(i) < k * h) )
      fun(i) = 4.0 * ( k * h  - x(i) ) * ( x(i) - (k-1)*h ) / h^2;
    else
      fun(i) = 0.0;
    end;
  end;
  
  

