[nder,xx,yfit,yp,n,x,c,work,ierr]=dpolvl(nder,xx,yfit,yp,n,x,c,work,ierr);

function [nder,xx,yfit,yp,n,x,c,work,ierr]=dpolvl(nder,xx,yfit,yp,n,x,c,work,ierr);
%***BEGIN PROLOGUE  DPOLVL
%***PURPOSE  Calculate the value of a polynomial and its first NDER
%            derivatives where the polynomial was produced by a previous
%            call to DPLINT.
%***LIBRARY   SLATEC
%***CATEGORY  E3
%***TYPE      doubleprecision (POLYVL-S, DPOLVL-D)
%***KEYWORDS  POLYNOMIAL EVALUATION
%***AUTHOR  Huddleston, R. E., (SNLL)
%***DESCRIPTION
%
%     Abstract -
%        subroutine DPOLVL calculates the value of the polynomial and
%     its first NDER derivatives where the polynomial was produced by
%     a previous call to DPLINT.
%        The variable N and the arrays X and C must not be altered
%     between the call to DPLINT and the call to DPOLVL.
%
%     ******  Dimensioning Information *******
%
%     YP   must be dimensioned by at least NDER
%     X    must be dimensioned by at least N (see the abstract )
%     C    must be dimensioned by at least N (see the abstract )
%     WORK must be dimensioned by at least 2*N if NDER is .GT. 0.
%
%     *** Note ***
%       If NDER=0, neither YP nor WORK need to be dimensioned variables.
%       If NDER=1, YP does not need to be a dimensioned variable.
%
%
%     *****  Input parameters
%       ***  All TYPE REAL variables are doubleprecision ***
%
%     NDER - the number of derivatives to be evaluated
%
%     XX   - the argument at which the polynomial and its derivatives
%            are to be evaluated.
%
%     N    - *****
%            *       N, X, and C must not be altered between the call
%     X    - *       to DPLINT and the call to DPOLVL.
%     C    - *****
%
%
%     *****  Output Parameters
%       ***  All TYPE REAL variables are doubleprecision ***
%
%     YFIT - the value of the polynomial at XX
%
%     YP   - the derivatives of the polynomial at XX.  The derivative of
%            order J at XX is stored in  YP(J) , J = 1,...,NDER.
%
%     IERR - Output error flag with the following possible values.
%          = 1  indicates normal execution
%
%     ***** Storage Parameters
%
%     WORK  = this is an array to provide internal working storage for
%             DPOLVL.  It must be dimensioned by at least 2*N if NDER is
%             .GT. 0.  If NDER=0, WORK does not need to be a dimensioned
%             variable.
%
%***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  (NONE)
%***REVISION HISTORY  (YYMMDD)
%   740601  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   891006  Cosmetic changes to prologue.  (WRB)
%   891006  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   920501  Reformatted the REFERENCES section.  (WRB)
%***end PROLOGUE  DPOLVL
persistent fac i im1 izero k km1 km1pi km2pn km2pni m mm ndr nmkp1 npkm1 pione pitwo pone ptwo xk ; 

if isempty(i), i=0; end;
if isempty(im1), im1=0; end;
if isempty(izero), izero=0; end;
if isempty(k), k=0; end;
if isempty(km1), km1=0; end;
if isempty(km1pi), km1pi=0; end;
if isempty(km2pn), km2pn=0; end;
if isempty(km2pni), km2pni=0; end;
if isempty(m), m=0; end;
if isempty(mm), mm=0; end;
if isempty(ndr), ndr=0; end;
if isempty(nmkp1), nmkp1=0; end;
if isempty(npkm1), npkm1=0; end;
c_shape=size(c);c=reshape(c,1,[]);
if isempty(fac), fac=0; end;
if isempty(pione), pione=0; end;
if isempty(pitwo), pitwo=0; end;
if isempty(pone), pone=0; end;
if isempty(ptwo), ptwo=0; end;
x_shape=size(x);x=reshape(x,1,[]);
if isempty(xk), xk=0; end;
yp_shape=size(yp);yp=reshape(yp,1,[]);
work_shape=size(work);work=reshape(work,1,[]);
%***FIRST EXECUTABLE STATEMENT  DPOLVL
ierr = 1;
if( nder<=0 )
%
%     *****   CODING FOR THE CASE NDER = 0
%
pione = 1.0d0;
pone = c(1);
yfit = pone;
if( n==1 )
c_shape=zeros(c_shape);c_shape(:)=c(1:numel(c_shape));c=c_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
yp_shape=zeros(yp_shape);yp_shape(:)=yp(1:numel(yp_shape));yp=yp_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
end;
for k = 2 : n;
pitwo =(xx-x(k-1)).*pione;
pione = pitwo;
ptwo = pone + pitwo.*c(k);
pone = ptwo;
end; k = fix(n+1);
yfit = ptwo;
c_shape=zeros(c_shape);c_shape(:)=c(1:numel(c_shape));c=c_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
yp_shape=zeros(yp_shape);yp_shape(:)=yp(1:numel(yp_shape));yp=yp_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
%
%     *****   end OF NDER = 0 CASE
%
elseif( n>1 ) ;
%
%     *****  end OF THE CASE  N = 1 AND  NDER .GT. 0
%
if( nder<n )
izero = 0;
ndr = fix(nder);
else;
%
%     *****  SET FLAGS FOR NUMBER OF DERIVATIVES AND FOR DERIVATIVES
%            IN EXCESS OF THE DEGREE (N-1) OF THE POLYNOMIAL.
%
izero = 1;
ndr = fix(n - 1);
end;
m = fix(ndr + 1);
mm = fix(m);
%
%     *****  START OF THE CASE NDER .GT. 0  AND N .GT. 1
%     *****  THE POLYNOMIAL AND ITS DERIVATIVES WILL BE EVALUATED AT XX
%
for k = 1 : ndr;
yp(k) = c(k+1);
end; k = fix(ndr+1);
%
%     *****  THE FOLLOWING SECTION OF CODE IS EASIER TO READ IF ONE
%            BREAKS WORK INTO TWO ARRAYS W AND V. THE CODE WOULD THEN
%            READ
%                W(1) = 1.
%                PONE = C(1)
%               *DO   K = 2,N
%               *   V(K-1) =  XX - X(K-1)
%               *   W(K)   =  V(K-1)*W(K-1)
%               *   PTWO   =  PONE + W(K)*C(K)
%               *   PONE   =  PWO
%
%               YFIT = PTWO
%
work(1) = 1.0d0;
pone = c(1);
for k = 2 : n;
km1 = fix(k - 1);
npkm1 = fix(n + k - 1);
work(npkm1) = xx - x(km1);
work(k) = work(npkm1).*work(km1);
ptwo = pone + work(k).*c(k);
pone = ptwo;
end; k = fix(n+1);
yfit = ptwo;
%
%     ** AT THIS POINT THE POLYNOMIAL HAS BEEN EVALUATED AND INFORMATION
%        FOR THE DERIVATIVE EVALUATIONS HAVE BEEN STORED IN THE ARRAY
%        WORK
if( n~=2 )
if( m==n )
mm = fix(ndr);
end;
%
%     ***** EVALUATE THE DERIVATIVES AT XX
%
%                  ******  DO K=2,MM   (FOR MOST CASES, MM = NDER + 1)
%                  *  ******  DO I=2,N-K+1
%                  *  *       W(I) = V(K-2+I)*W(I-1) + W(I)
%                  *  *       YP(K-1) = YP(K-1) + W(I)*C(K-1+I)
%                  ******  CONTINUE
%
for k = 2 : mm;
nmkp1 = fix(n - k + 1);
km1 = fix(k - 1);
km2pn = fix(k - 2 + n);
for i = 2 : nmkp1;
km2pni = fix(km2pn + i);
im1 = fix(i - 1);
km1pi = fix(km1 + i);
work(i) = work(km2pni).*work(im1) + work(i);
yp(km1) = yp(km1) + work(i).*c(km1pi);
end; i = fix(nmkp1+1);
end; k = fix(mm+1);
if( ndr~=1 )
fac = 1.0d0;
for k = 2 : ndr;
xk = k;
fac = xk.*fac;
yp(k) = fac.*yp(k);
end; k = fix(ndr+1);
end;
end;
%
%     ***** end OF DERIVATIVE EVALUATIONS
%
if( izero==0 )
c_shape=zeros(c_shape);c_shape(:)=c(1:numel(c_shape));c=c_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
yp_shape=zeros(yp_shape);yp_shape(:)=yp(1:numel(yp_shape));yp=yp_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
end;
%
%     *****  SET EXCESS DERIVATIVES TO ZERO.
%
for k = n : nder;
yp(k) = 0.0d0;
end; k = fix(nder+1);
else;
yfit = c(1);
%
%     *****  CODING FOR THE CASE  N=1 AND NDER .GT. 0
%
for k = 1 : nder;
yp(k) = 0.0d0;
end; k = fix(nder+1);
c_shape=zeros(c_shape);c_shape(:)=c(1:numel(c_shape));c=c_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
yp_shape=zeros(yp_shape);yp_shape(:)=yp(1:numel(yp_shape));yp=yp_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
end;
c_shape=zeros(c_shape);c_shape(:)=c(1:numel(c_shape));c=c_shape;
x_shape=zeros(x_shape);x_shape(:)=x(1:numel(x_shape));x=x_shape;
yp_shape=zeros(yp_shape);yp_shape(:)=yp(1:numel(yp_shape));yp=yp_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
end
%DECK DPOPT