function [] = plot_cubic (nel,x,Y,opt)
%--------------------------------------------------------
%  Plot C1 cubic between two nodal values of Y, & Y,x
%  nel = number of (two noded) elements
%  x   = mesh coordinates
%  Y   = computed waveform values at nodes
%  opt = plot type 1-Y 2-Y,x 3-Y & Y^2

n_gap = 10;  % number of intervals on curve for element
x_p = zeros (n_gap+1,1);   % local x value
Y_p = zeros (n_gap+1,1);   % local Y value
D_p = zeros (n_gap+1,1);   % local Y,x value
G_p = zeros (n_gap+1,2);   % local position interpolation
F_p = zeros (n_gap+1,4);   % local function interpolation
S_p = zeros (n_gap+1,4);   % local slope    interpolation

%  fill interpolation to use with each element (once)
for j = 1:n_gap+1
  r = 2*(j - 1)/n_gap - 1;  % local coordinate on -1 to +1
  G_p (j, 1) = 0.5*(1 - r) ;
  G_p (j, 2) = 0.5*(1 + r) ;
  F_p (j, 1) = (2 - 3*r + r^3)/4;          
  F_p (j, 2) = (1 - r - r^2 + r^3)/8;   % times element length
  F_p (j, 3) = (2 + 3*r - r^3)/4;    
  F_p (j, 4) = (-1 - r + r^2 + r^3)/8;  % times element length
  S_p (j, 1) = (-3 + 3*r*r) * 0.5;      % divide by element length
  S_p (j, 2) = (-1 - 2*r + 3*r*r) * 0.25;
  S_p (j, 3) = ( 3 - 3*r*r) * 0.5d0;    % divide by element length
  S_p (j, 4) = (-1 + 2*r + 3*r*r) * 0.25;
end % for j interior graph points

x_e   = zeros (2,1);   % element nodal x value
Y_e   = zeros (4,1);   % element nodal Y value
index = zeros (4,1);   % element nodal Y,x value

for iel = 1:nel    % loop for the total number of elements
  x_e (1) = x (iel);        % left x
  x_e (2) = x (iel+1);      % right x
  leng = x_e (2) - x_e (1); % element length

%   get element unknowns
  index = feeldof1 (iel,2,2);  % element dofs numbers
  Y_e = Y (index);             % element nodal Y1 Y1,x Y2 Y2,x

%   interpolate inside element & plot
  x_p (:) = G_p (:, 1:2) * x_e (1:2);               % mid points
  Y_p (:) = F_p (:, 1:2:4) * Y_e (1:2:4) ...
          + F_p (:, 2:2:4) * Y_e (2:2:4) * leng;    % mid values
  D_p (:) = S_p (:, 1:2:4) * Y_e (1:2:4) / leng ...
          + S_p (:, 2:2:4) * Y_e (2:2:4);           % mid slopes

%   plot this element part
  if ( opt == 1 )                  % plot Y only
    plot(x_e, Y_e(1:2:4),'ko')     % nodal value
    plot(x_p, Y_p,'k-')            % interpolated value
  elseif ( opt == 2 )              % plot Y,x only
    plot(x_e, Y_e(2:2:4),'kd')     % nodal slope
    plot(x_p, D_p,'k-')            % interpolated slope
  else                             % plot Y and Y^2
    plot(x_e, Y_e(1:2:4),'ko')     % nodal value
    plot(x_p, Y_p,'k-')            % interpolated value
    plot(x_e, Y_e(1:2:4).^2,'ks')  % nodal value
    plot(x_p, Y_p.^2,'k--')        % interpolated value
  end % if plot type
end % for iel elements
% end function plot_cubic