time   calls  line
                  1 function pc = pwch(x,y,s,dx,divdif)
                  2 %PWCH Piecewise cubic Hermite interpolation
                  3 %
                  4 %   PC = PWCH(X,Y,S) returns the ppform of the piecewise cubic Hermite
                  5 %   interpolant f to the given data values Y and the given slopes S
                  6 %   at the given data sites X. X must be a vector, and S must be of the same
                  7 %   size as Y.
                  8 %   If Y is a vector, then it must be of the same length as X and will be
                  9 %   resized if necessary to be a row vector.
                 10 %   If Y is a matrix, then size(Y,2) must equal length(X).
                 11 %
                 12 %   With d equal to size(Y,1), the piecewise cubic Hermite interpolant f to
                 13 %   these data is the d-valued piecewise cubic function with breaks at the
                 14 %   data sites X that satisfies
                 15 %
                 16 %          f(X(j)) = Y(:,j),  Df(X(j)) = S(:,j),  j=1:length(X) .
                 17 %
                 18 %   PC = PWCH(X,Y,S,DX,DIVDIF) also asks for DX := diff(X) and DIVDIF := 
                 19 %   diff(Y,1,2)./DX .
                 20 %
                 21 %   Class support for inputs X,Y,S,DX,DIVDIF:
                 22 %      float: double, single
                 23 %
                 24 %   See also SPLINE, INTERP1, PCHIP, PPVAL, UNMKPP.
                 25 
                 26 %   cb 1997
                 27 %   Copyright 1984-2010 The MathWorks, Inc.
                 28 
            25   29 if nargin<4, dx = diff(x(:).'); end 
  0.03      25   30 d = size(y,1); dxd = repmat(dx,d,1);  
            25   31 if nargin<5, divdif = diff(y,1,2)./dxd; end 
                 32 
            25   33 n = numel(x); 
  0.03      25   34 dzzdx = (divdif-s(:,1:n-1))./dxd; dzdxdx = (s(:,2:n)-divdif)./dxd; 
            25   35 dnm1 = d*(n-1); 
  0.10      25   36 pc = mkpp(x,[reshape((dzdxdx-dzzdx)./dxd,dnm1,1) ... 
                 37              reshape(2*dzzdx-dzdxdx,dnm1,1) ...
                 38              reshape(s(:,1:n-1),dnm1,1) ...
            25   39              reshape(y(:,1:n-1),dnm1,1)],d);