[dpchdfresult,k,x,s,ierr]=dpchdf(k,x,s,ierr);

function [dpchdfresult,k,x,s,ierr]=dpchdf(k,x,s,ierr);
dpchdfresult=[];
persistent firstCall i j value zero ; if isempty(firstCall),firstCall=1;end; 

;
%***BEGIN PROLOGUE  DPCHDF
%***SUBSIDIARY
%***PURPOSE  Computes divided differences for DPCHCE and DPCHSP
%***LIBRARY   SLATEC (PCHIP)
%***TYPE      doubleprecision (PCHDF-S, DPCHDF-D)
%***AUTHOR  Fritsch, F. N., (LLNL)
%***DESCRIPTION
%
%          DPCHDF:   DPCHIP Finite Difference Formula
%
%     Uses a divided difference formulation to compute a K-point approx-
%     imation to the derivative at X(K) based on the data in X and S.
%
%     Called by  DPCHCE  and  DPCHSP  to compute 3- and 4-point boundary
%     derivative approximations.
%
% ----------------------------------------------------------------------
%
%     On input:
%        K      is the order of the desired derivative approximation.
%               K must be at least 3 (error return if not).
%        X      contains the K values of the independent variable.
%               X need not be ordered, but the values **MUST** be
%               distinct.  (Not checked here.)
%        S      contains the associated slope values:
%                  S(I) = (F(I+1)-F(I))/(X(I+1)-X(I)), I=1(1)K-1.
%               (Note that S need only be of length K-1.)
%
%     On return:
%        S      will be destroyed.
%        IERR   will be set to -1 if K.LT.2 .
%        DPCHDF  will be set to the desired derivative approximation if
%               IERR=0 or to zero if IERR=-1.
%
% ----------------------------------------------------------------------
%
%***SEE ALSO  DPCHCE, DPCHSP
%***REFERENCES  Carl de Boor, A Practical Guide to Splines, Springer-
%                 Verlag, New York, 1978, pp. 10-16.
%***ROUTINES CALLED  XERMSG
%***REVISION HISTORY  (YYMMDD)
%   820503  DATE WRITTEN
%   820805  Converted to SLATEC library version.
%   870707  Corrected XERROR calls for d.p. name(s).
%   870813  Minor cosmetic changes.
%   890206  Corrected XERROR calls.
%   890411  Added SAVE statements (Vers. 3.2).
%   890411  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
%   900328  Added TYPE section.  (WRB)
%   910408  Updated AUTHOR and DATE WRITTEN sections in prologue.  (WRB)
%   920429  Revised format and order of references.  (WRB,FNF)
%   930503  Improved purpose.  (FNF)
%***end PROLOGUE  DPCHDF
%
%**end
%
%  DECLARE ARGUMENTS.
%
%
%  DECLARE LOCAL VARIABLES.
%
if isempty(i), i=0; end;
if isempty(j), j=0; end;
if isempty(value), value=0; end;
if isempty(zero), zero=0; end;
if firstCall,   zero=[0.0d0];  end;
firstCall=0;
%
%  CHECK FOR LEGAL VALUE OF K.
%
%***FIRST EXECUTABLE STATEMENT  DPCHDF
if( k<3 )
%
%  ERROR RETURN.
%
%     K.LT.3 RETURN.
ierr = -1;
[dumvar1,dumvar2,dumvar3,ierr]=xermsg('SLATEC','DPCHDF','K LESS THAN THREE',ierr,1);
dpchdfresult = zero;
else;
%
%  COMPUTE COEFFICIENTS OF INTERPOLATING POLYNOMIAL.
%
for j = 2 : k - 1;
for i = 1 : k - j;
s(i) =(s(i+1)-s(i))./(x(i+j)-x(i));
end; i = fix(k - j+1);
end; j = fix(k - 1+1);
%
%  EVALUATE DERIVATIVE AT X(K).
%
value = s(1);
for i = 2 : k - 1;
value = s(i) + value.*(x(k)-x(i));
end; i = fix(k - 1+1);
%
%  NORMAL RETURN.
%
ierr = 0;
dpchdfresult = value;
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(2)), assignin('caller','FUntemp',x); evalin('caller',[inputname(2),'=FUntemp;']); end
if csnil&&~isempty(inputname(3)), assignin('caller','FUntemp',s); evalin('caller',[inputname(3),'=FUntemp;']); end
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',k); evalin('caller',[inputname(1),'=FUntemp;']); end
if csnil&&~isempty(inputname(4)), assignin('caller','FUntemp',ierr); evalin('caller',[inputname(4),'=FUntemp;']); end
return;
end;
%------------- LAST LINE OF DPCHDF FOLLOWS -----------------------------
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(2)), assignin('caller','FUntemp',x); evalin('caller',[inputname(2),'=FUntemp;']); end
if csnil&&~isempty(inputname(3)), assignin('caller','FUntemp',s); evalin('caller',[inputname(3),'=FUntemp;']); end
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',k); evalin('caller',[inputname(1),'=FUntemp;']); end
if csnil&&~isempty(inputname(4)), assignin('caller','FUntemp',ierr); evalin('caller',[inputname(4),'=FUntemp;']); end
end
%DECK DPCHFD