function pp = ppfit2(x, xf, yf, method)
  a = zeros (length (xf), length (x)); 
  for ii = 1:length (x)
    p = zeros (length (x), 1);
    p(ii) = 1;
    a(:, ii) = interp1 (x(:), p, xf(:), method, "extrap");
  endfor
  pp = interp1 (x(:), a\yf(:), method, "pp");
endfunction

%!demo
%! clf
%! x = [0 1 2 3];
%! xf = linspace(0, 3, 1000);
%! yf = xf.^3 + rand(size(xf))*2.5;
%! plot (xf, yf, "y-")
%! hold on
%! pp = ppfit2(x,xf,yf, "linear");
%! plot (xf, ppval (pp, xf), "r-")
%! pp = ppfit2(x,xf,yf, "cubic");
%! plot (xf, ppval (pp, xf), "g-")
%! pp = ppfit2(x,xf,yf, "spline");
%! plot (xf, ppval (pp, xf), "b-")