function [xx,yy] = pwinterp(x,y,p,pts)
%PWINTERP Evaluate piecewise Lagrange interpolating polynomial
%
%	Given k*p+1 points (length(x) must be of the form k*p+1
%	for a positive integer k),
%
%	  x(1)<x(2)<...<x(k*p+1),
%
%	k*p+1 values y(1), y(2),...,y(k*p+1), an integer p>0,
%       and an integer pts>1, [xx,yy] = PWINTERP(x,y,p,pts) evaluates
%       the piecewise Lagrange interpolant of degree p on the k
%       intervals  (x(1),x(p+1)), (x(p+1),x(2*p+1)), ...
%	The vector xx gives the points where the interpolant is
%	evaluated.  These are the pts+1 evenly spaced points in
%	each subinterval, leading to k*pts+1 points altogether.
%	The vector yy gives the value of the interpolant at the
%	points in xx.
n=length(x);
if p < 1 | pts < 1 | rem(n-1,p) ~= 0 | n == 1 | n ~= length(y) 
  disp('error in input')
  return
end
nints = (n-1)/p ;
xbreak=x(1:p:n);
xx=zeros(nints*pts+1,1);
yy=xx;
j=1;
for int = 1:nints
  temp=linspace(xbreak(int),xbreak(int+1),pts+1);
  temp=temp(1:pts);
  xx(j:j+pts-1)=temp;
  yy(j:j+pts-1)=ninterp(x((int-1)*p+1:int*p+1),y((int-1)*p+1:int*p+1),temp);
  j=j+pts;
end
xx(j)=x(n);
yy(j)=y(n);

  
 
