function values = pcubic ( pc , XBAR )

%*****************************************************************************80
%
%MATLAB version of the SUBROUTINE PCUBIC on page 285 of conte/deboor
% Returns the value at XBAR  of the piecewise cubic function  f  described
% in the structure  pc .
% If  XBAR  is a sequence or even a higher-dimensional array,  values  will be
% of the same size as  XBAR, and will contain the value of  f  at 
% the corresponding element of  XBAR . 
% E.g., if XBAR == [1 2;3 4], then values == [f(1) f(2); f(3) f(4)].
%
%  Reference:
%
%    Samuel Conte, Carl de Boor,
%    Elementary Numerical Analysis,
%    Third Edition,
%    SIAM, 2017,
%    ISBN: 978-1-611975-19-2.
%
   XX = XBAR(:);  % enables the simultaneous evaluation of  f  at all the
                  % points in  XBAR  
   values = XBAR; % makes certain that the size of  values  matches  XBAR                 
XI = pc.breaks; nxi = length(XI);
   % locate  XX  within the break sequence  XI 
[~, index] = histc(XX,[-inf, XI, inf]); ii = max(2,min(nxi,index))-1;
   % this ensures that  XI(ii(j)) <= XX(j) < XI(ii(j)+1), except that,
   % for XX(j)<XI(2), ii(j)==1, and 
   % for XI(end-1)<=XX(j), ii(j)==nxi-1
   % hence using the first polynomial piece on (-inf .. XI(2)) and
   % using the last polynomial piece on [XI(nxi-1) .. inf) 
c = pc.coefs;
DX = XX.' - XI(ii); 
values(:) = c(1,ii) + DX.*(c(2,ii) + DX.*(c(3,ii) + DX.*c(4,ii)));
  return
end

