[maxord,mint,iswflg,el,tq]=ddcst(maxord,mint,iswflg,el,tq);

function [maxord,mint,iswflg,el,tq]=ddcst(maxord,mint,iswflg,el,tq);
%***BEGIN PROLOGUE  DDCST
%***SUBSIDIARY
%***PURPOSE  DDCST sets coefficients used by the core integrator DDSTP.
%***LIBRARY   SLATEC (SDRIVE)
%***TYPE      doubleprecision (SDCST-S, DDCST-D, CDCST-C)
%***AUTHOR  Kahaner, D. K., (NIST)
%             National Institute of Standards and Technology
%             Gaithersburg, MD  20899
%           Sutherland, C. D., (LANL)
%             Mail Stop D466
%             Los Alamos National Laboratory
%             Los Alamos, NM  87545
%***DESCRIPTION
%
%  DDCST is called by DDNTL.  The array EL determines the basic method.
%  The array TQ is involved in adjusting the step size in relation
%  to truncation error.  EL and TQ depend upon MINT, and are calculated
%  for orders 1 to MAXORD(.LE. 12).  For each order NQ, the coefficients
%  EL are calculated from the generating polynomial:
%    L(T) = EL(1,NQ) + EL(2,NQ)*T + ... + EL(NQ+1,NQ)*T**NQ.
%  For the implicit Adams methods, L(T) is given by
%    dL/dT = (1+T)*(2+T)* ... *(NQ-1+T)/K,   L(-1) = 0,
%    where      K = factorial(NQ-1).
%  For the Gear methods,
%    L(T) = (1+T)*(2+T)* ... *(NQ+T)/K,
%    where      K = factorial(NQ)*(1 + 1/2 + ... + 1/NQ).
%  For each order NQ, there are three components of TQ.
%
%***ROUTINES CALLED  (NONE)
%***REVISION HISTORY  (YYMMDD)
%   790601  DATE WRITTEN
%   900329  Initial submission to SLATEC.
%***end PROLOGUE  DDCST
persistent factrl gamma i j mxrd summlv ; 

el_orig=el;el_shape=[13,12];el=reshape([el_orig(1:min(prod(el_shape),numel(el_orig))),zeros(1,max(0,prod(el_shape)-numel(el_orig)))],el_shape);
if isempty(factrl), factrl=zeros(1,12); end;
if isempty(gamma), gamma=zeros(1,14); end;
if isempty(summlv), summlv=0; end;
tq_orig=tq;tq_shape=[3,12];tq=reshape([tq_orig(1:min(prod(tq_shape),numel(tq_orig))),zeros(1,max(0,prod(tq_shape)-numel(tq_orig)))],tq_shape);
if isempty(i), i=0; end;
if isempty(j), j=0; end;
if isempty(mxrd), mxrd=0; end;
%***FIRST EXECUTABLE STATEMENT  DDCST
factrl(1) = 1.0d0;
for i = 2 : maxord;
factrl(i) = i.*factrl(i-1);
end; i = fix(maxord+1);
%                                             Compute Adams coefficients
if( mint==1 )
gamma(1) = 1.0d0;
for i = 1 : maxord + 1;
summlv = 0.0d0;
for j = 1 : i;
summlv = summlv - gamma(j)./(i-j+2);
end; j = fix(i+1);
gamma(i+1) = summlv;
end; i = fix(maxord + 1+1);
el(1,1) = 1.0d0;
el(2,1) = 1.0d0;
el(2,2) = 1.0d0;
el(3,2) = 1.0d0;
for j = 3 : maxord;
el(2,j) = factrl(j-1);
for i = 3 : j;
el(i,j) =(j-1).*el(i,j-1) + el(i-1,j-1);
end; i = fix(j+1);
el(j+1,j) = 1.0d0;
end; j = fix(maxord+1);
for j = 2 : maxord;
el(1,j) = el(1,j-1) + gamma(j);
el(2,j) = 1.0d0;
for i = 3 : j + 1;
el(i,j) = el(i,j)./((i-1).*factrl(j-1));
end; i = fix(j + 1+1);
end; j = fix(maxord+1);
for j = 1 : maxord;
tq(1,j) = -1.0d0./(factrl(j).*gamma(j));
tq(2,j) = -1.0d0./gamma(j+1);
tq(3,j) = -1.0d0./gamma(j+2);
end; j = fix(maxord+1);
%                                              Compute Gear coefficients
elseif( mint==2 ) ;
el(1,1) = 1.0d0;
el(2,1) = 1.0d0;
for j = 2 : maxord;
el(1,j) = factrl(j);
for i = 2 : j;
el(i,j) = j.*el(i,j-1) + el(i-1,j-1);
end; i = fix(j+1);
el(j+1,j) = 1.0d0;
end; j = fix(maxord+1);
summlv = 1.0d0;
for j = 2 : maxord;
summlv = summlv + 1.0d0./j;
for i = 1 : j + 1;
el(i,j) = el(i,j)./(factrl(j).*summlv);
end; i = fix(j + 1+1);
end; j = fix(maxord+1);
for j = 1 : maxord;
if( j>1 )
tq(1,j) = 1.0d0./factrl(j-1);
end;
tq(2,j) =(j+1)./el(1,j);
tq(3,j) =(j+2)./el(1,j);
end; j = fix(maxord+1);
end;
%                          Compute constants used in the stiffness test.
%                          These are the ratio of TQ(2,NQ) for the Gear
%                          methods to those for the Adams methods.
if( iswflg==3 )
mxrd = fix(min(maxord,5));
if( mint==2 )
gamma(1) = 1.0d0;
for i = 1 : mxrd;
summlv = 0.0d0;
for j = 1 : i;
summlv = summlv - gamma(j)./(i-j+2);
end; j = fix(i+1);
gamma(i+1) = summlv;
end; i = fix(mxrd+1);
end;
summlv = 1.0d0;
for i = 2 : mxrd;
summlv = summlv + 1.0d0./i;
el(1+i,1) = -(i+1).*summlv.*gamma(i+1);
end; i = fix(mxrd+1);
end;
el_orig(1:prod(el_shape))=el;el=el_orig;
tq_orig(1:prod(tq_shape))=tq;tq=tq_orig;
end
%DECK DDEABM