function [xx,yy] = pwhermite(x,y,yp,pts)
%PWHERMITE Evaluate piecewise Hermite cubic interpolating polynomial
%
%	Given a vector of n points x (assumed to be distinct and increasing),
%	a vector of values y, and a vector of derivative values yp,
%       and a positive integer pts, [xx,yy]=PWHERMITE(x,y,yp,pts) returns
%       points xx and values yy, where xx contains pts+1 evenly spaced points
%       in each of the subintervals (x(i),x(i+1)) (so length(xx)=(n-1)*pts+1),
%       and yy contains the values of the piecewise Hermite cubic interpolant
%       and the points in xx.
n=length(x);
if n <2 | n ~= length(y) | n ~= length(yp)
  disp('error in input')
  return
end
xx=zeros((length(x)-1)*pts+1,1);
yy=xx;
j=1;
for i = 1:length(x)-1
  temp=linspace(x(i),x(i+1),pts+1);
  temp=temp(1:pts);
  xx(j:j+pts-1)=temp;
  yy(j:j+pts-1)=hermite3(x(i),x(i+1),y(i),y(i+1),yp(i),yp(i+1),temp);
  j=j+pts;
end
xx(j)=x(n);
yy(j)=y(n);

  
 
