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

function [meth,elco,tesco]=cfod(meth,elco,tesco);
%***BEGIN PROLOGUE  CFOD
%***SUBSIDIARY
%***PURPOSE  Subsidiary to DEBDF
%***LIBRARY   SLATEC
%***TYPE      SINGLE PRECISION (CFOD-S, DCFOD-D)
%***AUTHOR  (UNKNOWN)
%***DESCRIPTION
%
%   CFOD defines coefficients needed in the integrator package DEBDF
%
%***SEE ALSO  DEBDF
%***ROUTINES CALLED  (NONE)
%***REVISION HISTORY  (YYMMDD)
%   800901  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900328  Added TYPE section.  (WRB)
%***end PROLOGUE  CFOD
%
%
%LLL. OPTIMIZE
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(rqfac), rqfac=0; end;
if isempty(rq1fac), rq1fac=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);
%-----------------------------------------------------------------------
% CFOD  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.)
% CFOD  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  CFOD
if( meth==2 )
%
pc(1) = 1.0e0;
rq1fac = 1.0e0;
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.0e0;
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.0e0;
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.0e0;
elco(2,1) = 1.0e0;
tesco(1,1) = 0.0e0;
tesco(2,1) = 2.0e0;
tesco(1,2) = 1.0e0;
tesco(3,12) = 0.0e0;
pc(1) = 1.0e0;
rqfac = 1.0e0;
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.0e0;
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.0e0;
tsign = 1.0e0;
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.0e0;
for i = 2 : nq;
elco(i+1,nq) = rq1fac.*pc(i)./i;
end; i = fix(nq+1);
agamq = rqfac.*xpin;
ragq = 1.0e0./agamq;
tesco(2,nq) = ragq;
if( nq<12 )
tesco(1,nqp1) = ragq.*rqfac./nqp1;
end;
tesco(3,nqm1) = ragq;
end; nq = fix(12+1);
elco_orig(1:prod(elco_shape))=elco;elco=elco_orig;
tesco_orig(1:prod(tesco_shape))=tesco;tesco=tesco_orig;
return;
end;
%----------------------- end OF SUBROUTINE CFOD  -----------------------
elco_orig(1:prod(elco_shape))=elco;elco=elco_orig;
tesco_orig(1:prod(tesco_shape))=tesco;tesco=tesco_orig;
end
%DECK CGAMMA