[t,k,yh,nyh,dky,iflag]=dintyd(t,k,yh,nyh,dky,iflag);

function [t,k,yh,nyh,dky,iflag]=dintyd(t,k,yh,nyh,dky,iflag);
%***BEGIN PROLOGUE  DINTYD
%***SUBSIDIARY
%***PURPOSE  Subsidiary to DDEBDF
%***LIBRARY   SLATEC
%***TYPE      doubleprecision (INTYD-S, DINTYD-D)
%***AUTHOR  Watts, H. A., (SNLA)
%***DESCRIPTION
%
%   DINTYD approximates the solution and derivatives at T by polynomial
%   interpolation. Must be used in conjunction with the integrator
%   package DDEBDF.
% ----------------------------------------------------------------------
% DINTYD computes interpolated values of the K-th derivative of the
% dependent variable vector Y, and stores it in DKY.
% This routine is called by DDEBDF with K = 0,1 and T = TOUT, but may
% also be called by the user for any K up to the current order.
% (see detailed instructions in LSODE usage documentation.)
% ----------------------------------------------------------------------
% The computed values in DKY are gotten by interpolation using the
% Nordsieck history array YH.  This array corresponds uniquely to a
% vector-valued polynomial of degree NQCUR or less, and DKY is set
% to the K-th derivative of this polynomial at T.
% The formula for DKY is..
%              Q
%  DKY(I)  =  Sum  C(J,K) * (T - TN)**(J-K) * H**(-J) * YH(I,J+1)
%             J=K
% where  C(J,K) = J*(J-1)*...*(J-K+1), Q = NQCUR, TN = TCUR, H = HCUR.
% The quantities  NQ = NQCUR, L = NQ+1, N = NEQ, TN, and H are
% communicated by common.  The above sum is done in reverse order.
% IFLAG is returned negative if either K or T is out of bounds.
% ----------------------------------------------------------------------
%
%***SEE ALSO  DDEBDF
%***ROUTINES CALLED  (NONE)
%***COMMON BLOCKS    DDEBD1
%***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)
%   910722  Updated AUTHOR section.  (ALS)
%***end PROLOGUE  DINTYD
%
persistent c i ic j jb jb2 jj jj1 jp1 r s tp ; 

if isempty(i), i=0; end;
if isempty(ic), ic=0; end;
global ddebd1_12; if isempty(ddebd1_12), ddebd1_12=0; end;
global ddebd1_10; if isempty(ddebd1_10), ddebd1_10=zeros(1,14); end;
global ddebd1_11; if isempty(ddebd1_11), ddebd1_11=zeros(1,6); end;
if isempty(j), j=0; end;
if isempty(jb), jb=0; end;
if isempty(jb2), jb2=0; end;
if isempty(jj), jj=0; end;
if isempty(jj1), jj1=0; end;
if isempty(jp1), jp1=0; end;
global ddebd1_13; if isempty(ddebd1_13), ddebd1_13=0; end;
global ddebd1_14; if isempty(ddebd1_14), ddebd1_14=0; end;
global ddebd1_15; if isempty(ddebd1_15), ddebd1_15=0; end;
global ddebd1_18; if isempty(ddebd1_18), ddebd1_18=0; end;
global ddebd1_16; if isempty(ddebd1_16), ddebd1_16=0; end;
global ddebd1_17; if isempty(ddebd1_17), ddebd1_17=0; end;
global ddebd1_19; if isempty(ddebd1_19), ddebd1_19=0; end;
global ddebd1_22; if isempty(ddebd1_22), ddebd1_22=0; end;
global ddebd1_23; if isempty(ddebd1_23), ddebd1_23=0; end;
global ddebd1_20; if isempty(ddebd1_20), ddebd1_20=0; end;
global ddebd1_24; if isempty(ddebd1_24), ddebd1_24=0; end;
global ddebd1_21; if isempty(ddebd1_21), ddebd1_21=0; end;
if isempty(c), c=0; end;
global ddebd1_3; if isempty(ddebd1_3), ddebd1_3=0; end;
global ddebd1_4; if isempty(ddebd1_4), ddebd1_4=0; end;
global ddebd1_5; if isempty(ddebd1_5), ddebd1_5=0; end;
global ddebd1_6; if isempty(ddebd1_6), ddebd1_6=0; end;
global ddebd1_7; if isempty(ddebd1_7), ddebd1_7=0; end;
if isempty(r), r=0; end;
global ddebd1_1; if isempty(ddebd1_1), ddebd1_1=0; end;
global ddebd1_2; if isempty(ddebd1_2), ddebd1_2=zeros(1,210); end;
if isempty(s), s=0; end;
global ddebd1_8; if isempty(ddebd1_8), ddebd1_8=0; end;
if isempty(tp), tp=0; end;
global ddebd1_9; if isempty(ddebd1_9), ddebd1_9=0; end;
yh_shape=size(yh);yh=reshape([yh(:).',zeros(1,ceil(numel(yh)./prod([nyh])).*prod([nyh])-numel(yh))],nyh,[]);
dky_shape=size(dky);dky=reshape(dky,1,[]);
% common :: ;
%% common /ddebd1/ rownd , rowns(210) , el0 , h , hmin , hmxi , hu ,tn , uround , iownd(14) , iowns(6) , ier ,jstart , kflag , l , meth , miter , maxord , n ,nq , nst , nfe , nje , nqu;
%% common /ddebd1/ ddebd1_1 , ddebd1_2(210) , ddebd1_3 , ddebd1_4 , ddebd1_5 , ddebd1_6 , ddebd1_7 ,ddebd1_8 , ddebd1_9 , ddebd1_10(14) , ddebd1_11(6) , ddebd1_12 ,ddebd1_13 , ddebd1_14 , ddebd1_15 , ddebd1_16 , ddebd1_17 , ddebd1_18 , ddebd1_19 ,ddebd1_20 , ddebd1_21 , ddebd1_22 , ddebd1_23 , ddebd1_24;
%
%     BEGIN BLOCK PERMITTING ...EXITS TO 130
%***FIRST EXECUTABLE STATEMENT  DINTYD
iflag = 0;
if( k<0 || k>ddebd1_20 )
%
iflag = -1;
else;
tp = ddebd1_8 - ddebd1_7.*(1.0d0+100.0d0.*ddebd1_9);
if((t-tp).*(t-ddebd1_8)<=0.0d0 )
%
s =(t-ddebd1_8)./ddebd1_4;
ic = 1;
if( k~=0 )
jj1 = fix(ddebd1_15 - k);
for jj = jj1 : ddebd1_20;
ic = fix(ic.*jj);
end; jj = fix(ddebd1_20+1);
end;
c = ic;
for i = 1 : ddebd1_19;
dky(i) = c.*yh(i,ddebd1_15);
end; i = fix(ddebd1_19+1);
if( k~=ddebd1_20 )
jb2 = fix(ddebd1_20 - k);
for jb = 1 : jb2;
j = fix(ddebd1_20 - jb);
jp1 = fix(j + 1);
ic = 1;
if( k~=0 )
jj1 = fix(jp1 - k);
for jj = jj1 : j;
ic = fix(ic.*jj);
end; jj = fix(j+1);
end;
c = ic;
for i = 1 : ddebd1_19;
dky(i) = c.*yh(i,jp1) + s.*dky(i);
end; i = fix(ddebd1_19+1);
end; jb = fix(jb2+1);
%     .........EXIT
if( k==0 )
yh_shape=zeros(yh_shape);yh_shape(:)=yh(1:numel(yh_shape));yh=yh_shape;
dky_shape=zeros(dky_shape);dky_shape(:)=dky(1:numel(dky_shape));dky=dky_shape;
return;
end;
end;
r = ddebd1_4.^(-k);
for i = 1 : ddebd1_19;
dky(i) = r.*dky(i);
end; i = fix(ddebd1_19+1);
else;
%     .........EXIT
iflag = -2;
end;
end;
%     ----------------------- end OF SUBROUTINE DINTYD
%     -----------------------
yh_shape=zeros(yh_shape);yh_shape(:)=yh(1:numel(yh_shape));yh=yh_shape;
dky_shape=zeros(dky_shape);dky_shape(:)=dky(1:numel(dky_shape));dky=dky_shape;
end
%DECK DIR