[l,nder,x,yfit,yp,a]=dp1vlu(l,nder,x,yfit,yp,a);

function [l,nder,x,yfit,yp,a]=dp1vlu(l,nder,x,yfit,yp,a);
%***BEGIN PROLOGUE  DP1VLU
%***PURPOSE  Use the coefficients generated by DPOLFT to evaluate the
%            polynomial fit of degree L, along with the first NDER of
%            its derivatives, at a specified point.
%***LIBRARY   SLATEC
%***CATEGORY  K6
%***TYPE      doubleprecision (PVALUE-S, DP1VLU-D)
%***KEYWORDS  CURVE FITTING, LEAST SQUARES, POLYNOMIAL APPROXIMATION
%***AUTHOR  Shampine, L. F., (SNLA)
%           Davenport, S. M., (SNLA)
%***DESCRIPTION
%
%     Abstract
%
%     The subroutine  DP1VLU  uses the coefficients generated by  DPOLFT
%     to evaluate the polynomial fit of degree  L , along with the first
%     NDER  of its derivatives, at a specified point.  Computationally
%     stable recurrence relations are used to perform this task.
%
%     The parameters for  DP1VLU  are
%
%     Input -- ALL TYPE REAL variables are doubleprecision
%         L -      the degree of polynomial to be evaluated.  L  may be
%                  any non-negative integer which is less than or equal
%                  to  NDEG , the highest degree polynomial provided
%                  by  DPOLFT .
%         NDER -   the number of derivatives to be evaluated.  NDER
%                  may be 0 or any positive value.  If NDER is less
%                  than 0, it will be treated as 0.
%         X -      the argument at which the polynomial and its
%                  derivatives are to be evaluated.
%         A -      work and output array containing values from last
%                  call to  DPOLFT .
%
%     Output -- ALL TYPE REAL variables are doubleprecision
%         YFIT -   value of the fitting polynomial of degree  L  at  X
%         YP -     array containing the first through  NDER  derivatives
%                  of the polynomial of degree  L .  YP  must be
%                  dimensioned at least  NDER  in the calling program.
%
%***REFERENCES  L. F. Shampine, S. M. Davenport and R. E. Huddleston,
%                 Curve fitting by polynomials in one variable, Report
%                 SLA-74-0270, Sandia Laboratories, June 1974.
%***ROUTINES CALLED  XERMSG
%***REVISION HISTORY  (YYMMDD)
%   740601  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890911  Removed unnecessary intrinsics.  (WRB)
%   891006  Cosmetic changes to prologue.  (WRB)
%   891006  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)
%   900510  Convert XERRWV calls to XERMSG calls.  (RWC)
%   920501  Reformatted the REFERENCES section.  (WRB)
%***end PROLOGUE  DP1VLU
persistent cc dif i ic ilo in inp1 iup k1 k1i k2 k3 k3p1 k3pn k4 k4p1 k4pn kc lm1 lp1 maxord n ndo ndp1 nord val xern1 xern2 ; 

if isempty(i), i=0; end;
if isempty(ic), ic=0; end;
if isempty(ilo), ilo=0; end;
if isempty(in), in=0; end;
if isempty(inp1), inp1=0; end;
if isempty(iup), iup=0; end;
if isempty(k1), k1=0; end;
if isempty(k1i), k1i=0; end;
if isempty(k2), k2=0; end;
if isempty(k3), k3=0; end;
if isempty(k3p1), k3p1=0; end;
if isempty(k3pn), k3pn=0; end;
if isempty(k4), k4=0; end;
if isempty(k4p1), k4p1=0; end;
if isempty(k4pn), k4pn=0; end;
if isempty(kc), kc=0; end;
if isempty(lm1), lm1=0; end;
if isempty(lp1), lp1=0; end;
if isempty(maxord), maxord=0; end;
if isempty(n), n=0; end;
if isempty(ndo), ndo=0; end;
if isempty(ndp1), ndp1=0; end;
if isempty(nord), nord=0; end;
a_shape=size(a);a=reshape(a,1,[]);
if isempty(cc), cc=0; end;
if isempty(dif), dif=0; end;
if isempty(val), val=0; end;
yp_shape=size(yp);yp=reshape(yp,1,[]);
if isempty(xern1), xern1=repmat(' ',1,8); end;
if isempty(xern2), xern2=repmat(' ',1,8); end;
%***FIRST EXECUTABLE STATEMENT  DP1VLU
if( l<0 )
%
xermsg('SLATEC','DP1VLU',['INVALID INPUT PARAMETER.  ORDER OF POLYNOMIAL EVALUATION ','REQUESTED IS NEGATIVE.'],2,2);
else;
ndo = fix(max(nder,0));
ndo = fix(min(ndo,l));
maxord = fix(a(1) + 0.5d0);
k1 = fix(maxord + 1);
k2 = fix(k1 + maxord);
k3 = fix(k2 + maxord + 2);
nord = fix(a(k3) + 0.5d0);
if( l>nord )
%
xern1=sprintf(['%8i'], l);
xern2=sprintf(['%8i'], nord);
xermsg('SLATEC','DP1VLU',['THE ORDER OF POLYNOMIAL EVALUATION, L = ',[xern1,[' REQUESTED EXCEEDS THE HIGHEST ORDER FIT, NORD = ',[xern2,', COMPUTED BY DPOLFT -- EXECUTION TERMINATED.']]]],8,2);
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
yp_shape=zeros(yp_shape);yp_shape(:)=yp(1:numel(yp_shape));yp=yp_shape;
return;
else;
k4 = fix(k3 + l + 1);
if( nder>=1 )
for i = 1 : nder;
yp(i) = 0.0d0;
end; i = fix(nder+1);
end;
if( l>=2 )
%
% L IS GREATER THAN 1
%
ndp1 = fix(ndo + 1);
k3p1 = fix(k3 + 1);
k4p1 = fix(k4 + 1);
lp1 = fix(l + 1);
lm1 = fix(l - 1);
ilo = fix(k3 + 3);
iup = fix(k4 + ndp1);
for i = ilo : iup;
a(i) = 0.0d0;
end; i = fix(iup+1);
dif = x - a(lp1);
kc = fix(k2 + lp1);
a(k4p1) = a(kc);
a(k3p1) = a(kc-1) + dif.*a(k4p1);
a(k3+2) = a(k4p1);
%
% EVALUATE RECURRENCE RELATIONS FOR FUNCTION VALUE AND DERIVATIVES
%
for i = 1 : lm1;
in = fix(l - i);
inp1 = fix(in + 1);
k1i = fix(k1 + inp1);
ic = fix(k2 + in);
dif = x - a(inp1);
val = a(ic) + dif.*a(k3p1) - a(k1i).*a(k4p1);
if( ndo>0 )
for n = 1 : ndo;
k3pn = fix(k3p1 + n);
k4pn = fix(k4p1 + n);
yp(n) = dif.*a(k3pn) + n.*a(k3pn-1) - a(k1i).*a(k4pn);
end; n = fix(ndo+1);
%
% SAVE VALUES NEEDED FOR NEXT EVALUATION OF RECURRENCE RELATIONS
%
for n = 1 : ndo;
k3pn = fix(k3p1 + n);
k4pn = fix(k4p1 + n);
a(k4pn) = a(k3pn);
a(k3pn) = yp(n);
end; n = fix(ndo+1);
end;
a(k4p1) = a(k3p1);
a(k3p1) = val;
end; i = fix(lm1+1);
elseif( l==1 ) ;
%
% L IS 1
%
cc = a(k2+2);
val = a(k2+1) +(x-a(2)).*cc;
if( nder>=1 )
yp(1) = cc;
end;
else;
%
% L IS 0
%
val = a(k2+1);
end;
%
% NORMAL RETURN OR ABORT DUE TO ERROR
%
yfit = val;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
yp_shape=zeros(yp_shape);yp_shape(:)=yp(1:numel(yp_shape));yp=yp_shape;
return;
end;
end;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
yp_shape=zeros(yp_shape);yp_shape(:)=yp(1:numel(yp_shape));yp=yp_shape;
end
%DECK DPBCO