[meth,elco,tesco]=dcfod(meth,elco,tesco);

function [meth,elco,tesco]=dcfod(meth,elco,tesco);
%***BEGIN PROLOGUE  DCFOD
%***SUBSIDIARY
%***PURPOSE  Subsidiary to DDEBDF
%***LIBRARY   SLATEC
%***TYPE      doubleprecision (CFOD-S, DCFOD-D)
%***AUTHOR  (UNKNOWN)
%***DESCRIPTION
%
%   DCFOD defines coefficients needed in the integrator package DDEBDF
%
%***SEE ALSO  DDEBDF
%***ROUTINES CALLED  (NONE)
%***REVISION HISTORY  (YYMMDD)
%   820301  DATE WRITTEN
%   890911  Removed unnecessary intrinsics.  (WRB)
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900328  Added TYPE section.  (WRB)
%***end PROLOGUE  DCFOD
%
%
persistent agamq fnq fnqm1 i ib nq nqm1 nqp1 pc pint ragq rq1fac rqfac tsign xpin ; 

if isempty(i), i=0; end;
if isempty(ib), ib=0; end;
if isempty(nq), nq=0; end;
if isempty(nqm1), nqm1=0; end;
if isempty(nqp1), nqp1=0; end;
if isempty(agamq), agamq=0; end;
if isempty(fnq), fnq=0; end;
if isempty(fnqm1), fnqm1=0; end;
if isempty(pc), pc=zeros(1,12); end;
if isempty(pint), pint=0; end;
if isempty(ragq), ragq=0; end;
if isempty(rq1fac), rq1fac=0; end;
if isempty(rqfac), rqfac=0; end;
if isempty(tsign), tsign=0; end;
if isempty(xpin), xpin=0; end;
elco_orig=elco;elco_shape=[13,12];elco=reshape([elco_orig(1:min(prod(elco_shape),numel(elco_orig))),zeros(1,max(0,prod(elco_shape)-numel(elco_orig)))],elco_shape);
tesco_orig=tesco;tesco_shape=[3,12];tesco=reshape([tesco_orig(1:min(prod(tesco_shape),numel(tesco_orig))),zeros(1,max(0,prod(tesco_shape)-numel(tesco_orig)))],tesco_shape);
%     ------------------------------------------------------------------
%      DCFOD  IS CALLED BY THE INTEGRATOR ROUTINE TO SET COEFFICIENTS
%      NEEDED THERE.  THE COEFFICIENTS FOR THE CURRENT METHOD, AS
%      GIVEN BY THE VALUE OF METH, ARE SET FOR ALL ORDERS AND SAVED.
%      THE MAXIMUM ORDER ASSUMED HERE IS 12 IF METH = 1 AND 5 IF METH =
%      2.  (A SMALLER VALUE OF THE MAXIMUM ORDER IS ALSO ALLOWED.)
%      DCFOD  IS CALLED ONCE AT THE BEGINNING OF THE PROBLEM,
%      AND IS NOT CALLED AGAIN UNLESS AND UNTIL METH IS CHANGED.
%
%      THE ELCO ARRAY CONTAINS THE BASIC METHOD COEFFICIENTS.
%      THE COEFFICIENTS EL(I), 1 .LE. I .LE. NQ+1, FOR THE METHOD OF
%      ORDER NQ ARE STORED IN ELCO(I,NQ).  THEY ARE GIVEN BY A
%      GENERATING POLYNOMIAL, I.E.,
%          L(X) = EL(1) + EL(2)*X + ... + EL(NQ+1)*X**NQ.
%      FOR THE IMPLICIT ADAMS METHODS, L(X) IS GIVEN BY
%          DL/DX = (X+1)*(X+2)*...*(X+NQ-1)/FACTORIAL(NQ-1),    L(-1) =
%      0.  FOR THE BDF METHODS, L(X) IS GIVEN BY
%          L(X) = (X+1)*(X+2)* ... *(X+NQ)/K,
%      WHERE         K = FACTORIAL(NQ)*(1 + 1/2 + ... + 1/NQ).
%
%      THE TESCO ARRAY CONTAINS TEST CONSTANTS USED FOR THE
%      LOCAL ERROR TEST AND THE SELECTION OF STEP SIZE AND/OR ORDER.
%      AT ORDER NQ, TESCO(K,NQ) IS USED FOR THE SELECTION OF STEP
%      SIZE AT ORDER NQ - 1 IF K = 1, AT ORDER NQ IF K = 2, AND AT ORDER
%      NQ + 1 IF K = 3.
%     ------------------------------------------------------------------
%
%***FIRST EXECUTABLE STATEMENT  DCFOD
if( meth==2 )
%
pc(1) = 1.0d0;
rq1fac = 1.0d0;
for nq = 1 : 5;
%           ------------------------------------------------------------
%            THE PC ARRAY WILL CONTAIN THE COEFFICIENTS OF THE
%                POLYNOMIAL P(X) = (X+1)*(X+2)*...*(X+NQ).
%            INITIALLY, P(X) = 1.
%           ------------------------------------------------------------
fnq = nq;
nqp1 = fix(nq + 1);
%           FORM COEFFICIENTS OF P(X)*(X+NQ).
%           ------------------------------------
pc(nqp1) = 0.0d0;
for ib = 1 : nq;
i = fix(nq + 2 - ib);
pc(i) = pc(i-1) + fnq.*pc(i);
end; ib = fix(nq+1);
pc(1) = fnq.*pc(1);
%           STORE COEFFICIENTS IN ELCO AND TESCO.
%           --------------------------------
for i = 1 : nqp1;
elco(i,nq) = pc(i)./pc(2);
end; i = fix(nqp1+1);
elco(2,nq) = 1.0d0;
tesco(1,nq) = rq1fac;
tesco(2,nq) = nqp1./elco(1,nq);
tesco(3,nq) =(nq+2)./elco(1,nq);
rq1fac = rq1fac./fnq;
end; nq = fix(5+1);
else;
%
elco(1,1) = 1.0d0;
elco(2,1) = 1.0d0;
tesco(1,1) = 0.0d0;
tesco(2,1) = 2.0d0;
tesco(1,2) = 1.0d0;
tesco(3,12) = 0.0d0;
pc(1) = 1.0d0;
rqfac = 1.0d0;
for nq = 2 : 12;
%           ------------------------------------------------------------
%            THE PC ARRAY WILL CONTAIN THE COEFFICIENTS OF THE
%                POLYNOMIAL P(X) = (X+1)*(X+2)*...*(X+NQ-1).
%            INITIALLY, P(X) = 1.
%           ------------------------------------------------------------
rq1fac = rqfac;
rqfac = rqfac./nq;
nqm1 = fix(nq - 1);
fnqm1 = nqm1;
nqp1 = fix(nq + 1);
%           FORM COEFFICIENTS OF P(X)*(X+NQ-1).
%           ----------------------------------
pc(nq) = 0.0d0;
for ib = 1 : nqm1;
i = fix(nqp1 - ib);
pc(i) = pc(i-1) + fnqm1.*pc(i);
end; ib = fix(nqm1+1);
pc(1) = fnqm1.*pc(1);
%           COMPUTE INTEGRAL, -1 TO 0, OF P(X) AND X*P(X).
%           -----------------------
pint = pc(1);
xpin = pc(1)./2.0d0;
tsign = 1.0d0;
for i = 2 : nq;
tsign = -tsign;
pint = pint + tsign.*pc(i)./i;
xpin = xpin + tsign.*pc(i)./(i+1);
end; i = fix(nq+1);
%           STORE COEFFICIENTS IN ELCO AND TESCO.
%           --------------------------------
elco(1,nq) = pint.*rq1fac;
elco(2,nq) = 1.0d0;
for i = 2 : nq;
elco(i+1,nq) = rq1fac.*pc(i)./i;
end; i = fix(nq+1);
agamq = rqfac.*xpin;
ragq = 1.0d0./agamq;
tesco(2,nq) = ragq;
if( nq<12 )
tesco(1,nqp1) = ragq.*rqfac./nqp1;
end;
tesco(3,nqm1) = ragq;
end; nq = fix(12+1);
end;
%     ----------------------- end OF SUBROUTINE DCFOD
%     -----------------------
elco_orig(1:prod(elco_shape))=elco;elco=elco_orig;
tesco_orig(1:prod(tesco_shape))=tesco;tesco=tesco_orig;
end
%DECK DCG