function [v, vp, vpp] = pqhermite(x, yin, u)
%PQHERMITE  True piecewise quintic Hermite interpolation.
%  [v, vp, vpp] = pqhip(x,y,u) evaluates the true piecewise quintic
%  interpolant P(x), with P(x(j)) = y(1,j), P'(x(j)) = y(2,j), 
%  P''(x(j)) = y(3,j)
%  and returns v(k) = P(u(k)), vp(k) = P'(u(k)), and vpp(k) = P''(u(k)).
%
%  x = sort(x);  nodes  must be sorted 
%  RMC November 2012 modified from pqhip
 
% -------------------------------------------------------

   %x = sort(x);  nodes  must be sorted 
   h = diff(x);
   y(1:3,:) = yin(1:3,:);



%  Piecewise polynomial coefficients are just the data themselves,
%  because we use the Hermite interpolational basis.

   n = length(x);

%  Find subinterval indices k so that x(k) <= u < x(k+1)

   k = ones(size(u));
   for j = 2:n-1
      k(x(j) <= u) = j;
   end

%  Evaluate interpolant, using formulas generated by Maple
%  and modified by hand: * --> .*, ^ --> .^, and / --> ./

   s = u - x(k);
   h = x(k+1)-x(k);
   if any( h <= 0),
       error 'This routine requires strictly increasing real nodes.';
   end
   theta = s./h;

   v = -((theta - 1) .^ 3 .* (6 * theta .^ 2 + 1 + 3 * theta) .* y(1, k)) ...
       - (theta .* h .* (theta - 1) .^ 3 .* (1 + 3 * theta) .* y(2, k)) ...
       - (theta .^ 2 .* h .^ 2 .* (theta - 1) .^ 3 .* y(3, k)) / 0.2e1 ...
       + (theta .^ 3 .* (6 * theta .^ 2 - 15 * theta + 10) .* y(1, k + 1)) ...
       - (theta .^ 3 .* h .* (theta - 1) .* (3 * theta - 4) .* y(2, k + 1)) ...
       + (theta .^ 3 .* h .^ 2 .* (theta - 1) .^ 2 .* y(3, k + 1)) / 0.2e1;

    vp = -(30 * theta .^ 2 .* (theta - 1) .^ 2 .* y(1, k))./h ...
        - ((theta - 1) .^ 2 .* (5 * theta + 1) .* (3 * theta - 1) .* y(2, k)) ...
        - (theta .* h .* (theta - 1) .^ 2 .* (5 * theta - 2) .* y(3, k)) / 0.2e1 ...
        + (30 .* theta .^ 2 .* (theta - 1) .^ 2 .* y(1, k + 1))./h ...
        - (theta .^ 2 .* (3 * theta - 2) .* (5 * theta - 6) .* y(2, k + 1)) ...
        + (theta .^ 2 .* h .* (theta - 1) .* (5 * theta - 3) .* y(3, k + 1)) / 0.2e1;

    vpp = 60  * (theta - 1) .* theta .* (2 * theta - 1) .* (y(1,k+1)-y(1, k))./ h .^ 2 ...
        - 12 * (theta - 1) .* theta .* (5 * theta - 3) .* y(2, k)./h ...
        - (theta - 1) .* (10 * theta .^ 2 - 8 * theta + 1) .* y(3, k) ...
        - 12 * theta .* (5 * theta - 2) .* (theta - 1) .* y(2, k + 1)./h ...
        + theta .* (10 * theta .^ 2 - 12 * theta + 3) .* y(3, k + 1);
end