[dcsevlresult,x,cs,n]=dcsevl(x,cs,n);

function [dcsevlresult,x,cs,n]=dcsevl(x,cs,n);
dcsevlresult=[];
persistent b0 b1 b2 first firstCall i ni onepl twox ; if isempty(firstCall),firstCall=1;end; 

;
if isempty(i), i=0; end;
if isempty(ni), ni=0; end;
%***BEGIN PROLOGUE  DCSEVL
%***PURPOSE  Evaluate a Chebyshev series.
%***LIBRARY   SLATEC (FNLIB)
%***CATEGORY  C3A2
%***TYPE      doubleprecision (CSEVL-S, DCSEVL-D)
%***KEYWORDS  CHEBYSHEV SERIES, FNLIB, SPECIAL FUNCTIONS
%***AUTHOR  Fullerton, W., (LANL)
%***DESCRIPTION
%
%  Evaluate the N-term Chebyshev series CS at X.  Adapted from
%  a method presented in the paper by Broucke referenced below.
%
%       Input Arguments --
%  X    value at which the series is to be evaluated.
%  CS   array of N terms of a Chebyshev series.  In evaluating
%       CS, only half the first coefficient is summed.
%  N    number of terms in array CS.
%
%***REFERENCES  R. Broucke, Ten subroutines for the manipulation of
%                 Chebyshev series, Algorithm 446, Communications of
%                 the A.C.M. 16, (1973) pp. 254-256.
%               L. Fox and I. B. Parker, Chebyshev Polynomials in
%                 Numerical Analysis, Oxford University Press, 1968,
%                 page 56.
%***ROUTINES CALLED  D1MACH, XERMSG
%***REVISION HISTORY  (YYMMDD)
%   770401  DATE WRITTEN
%   890831  Modified array declarations.  (WRB)
%   890831  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)
%   900329  Prologued revised extensively and code rewritten to allow
%           X to be slightly outside interval (-1,+1).  (WRB)
%   920501  Reformatted the REFERENCES section.  (WRB)
%***end PROLOGUE  DCSEVL
if isempty(b0), b0=0; end;
if isempty(b1), b1=0; end;
if isempty(b2), b2=0; end;
cs_shape=size(cs);cs=reshape(cs,1,[]);
if isempty(onepl), onepl=0; end;
if isempty(twox), twox=0; end;
if isempty(first), first=false; end;
if firstCall,   first=[true];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  DCSEVL
if( first )
onepl = 1.0d0 + d1mach(4);
end;
first = false;
if( n<1 )
xermsg('SLATEC','DCSEVL','NUMBER OF TERMS .LE. 0',2,2);
end;
if( n>1000 )
xermsg('SLATEC','DCSEVL','NUMBER OF TERMS .GT. 1000',3,2);
end;
if( abs(x)>onepl )
xermsg('SLATEC','DCSEVL','X OUTSIDE THE INTERVAL (-1,+1)',1,1);
end;
%
b1 = 0.0d0;
b0 = 0.0d0;
twox = 2.0d0.*x;
for i = 1 : n;
b2 = b1;
b1 = b0;
ni = fix(n + 1 - i);
b0 = twox.*b1 - b2 + cs(ni);
end; i = fix(n+1);
%
dcsevlresult = 0.5d0.*(b0-b2);
%
cs_shape=zeros(cs_shape);cs_shape(:)=cs(1:numel(cs_shape));cs=cs_shape;
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',x); evalin('caller',[inputname(1),'=FUntemp;']); end
if csnil&&~isempty(inputname(3)), assignin('caller','FUntemp',n); evalin('caller',[inputname(3),'=FUntemp;']); end
if csnil&&~isempty(inputname(2)), assignin('caller','FUntemp',cs); evalin('caller',[inputname(2),'=FUntemp;']); end
end
%DECK DCV