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

function [t,k,yh,nyh,dky,iflag]=intyd(t,k,yh,nyh,dky,iflag);
%***BEGIN PROLOGUE  INTYD
%***SUBSIDIARY
%***PURPOSE  Subsidiary to DEBDF
%***LIBRARY   SLATEC
%***TYPE      SINGLE PRECISION (INTYD-S, DINTYD-D)
%***AUTHOR  Watts, H. A., (SNLA)
%***DESCRIPTION
%
%   INTYD approximates the solution and derivatives at T by polynomial
%   interpolation. Must be used in conjunction with the integrator
%   package DEBDF.
% ----------------------------------------------------------------------
% INTYD computes interpolated values of the K-th derivative of the
% dependent variable vector Y, and stores it in DKY.
% This routine is called by DEBDF 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  DEBDF
%***ROUTINES CALLED  (NONE)
%***COMMON BLOCKS    DEBDF1
%***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)
%   910722  Updated AUTHOR section.  (ALS)
%***end PROLOGUE  INTYD
%
%LLL. OPTIMIZE
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 debdf1_12; if isempty(debdf1_12), debdf1_12=0; end;
global debdf1_10; if isempty(debdf1_10), debdf1_10=zeros(1,14); end;
global debdf1_11; if isempty(debdf1_11), debdf1_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 debdf1_13; if isempty(debdf1_13), debdf1_13=0; end;
global debdf1_14; if isempty(debdf1_14), debdf1_14=0; end;
global debdf1_15; if isempty(debdf1_15), debdf1_15=0; end;
global debdf1_18; if isempty(debdf1_18), debdf1_18=0; end;
global debdf1_16; if isempty(debdf1_16), debdf1_16=0; end;
global debdf1_17; if isempty(debdf1_17), debdf1_17=0; end;
global debdf1_19; if isempty(debdf1_19), debdf1_19=0; end;
global debdf1_22; if isempty(debdf1_22), debdf1_22=0; end;
global debdf1_23; if isempty(debdf1_23), debdf1_23=0; end;
global debdf1_20; if isempty(debdf1_20), debdf1_20=0; end;
global debdf1_24; if isempty(debdf1_24), debdf1_24=0; end;
global debdf1_21; if isempty(debdf1_21), debdf1_21=0; end;
global debdf1_1; if isempty(debdf1_1), debdf1_1=0; end;
global debdf1_2; if isempty(debdf1_2), debdf1_2=zeros(1,210); end;
global debdf1_3; if isempty(debdf1_3), debdf1_3=0; end;
global debdf1_4; if isempty(debdf1_4), debdf1_4=0; end;
global debdf1_5; if isempty(debdf1_5), debdf1_5=0; end;
global debdf1_6; if isempty(debdf1_6), debdf1_6=0; end;
global debdf1_7; if isempty(debdf1_7), debdf1_7=0; end;
global debdf1_8; if isempty(debdf1_8), debdf1_8=0; end;
global debdf1_9; if isempty(debdf1_9), debdf1_9=0; end;
if isempty(c), c=0; end;
if isempty(r), r=0; end;
if isempty(s), s=0; end;
if isempty(tp), tp=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 /debdf1/ 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 /debdf1/ debdf1_1 , debdf1_2(210) , debdf1_3 , debdf1_4 , debdf1_5 , debdf1_6 , debdf1_7 ,debdf1_8 , debdf1_9 , debdf1_10(14) , debdf1_11(6) , debdf1_12 ,debdf1_13 , debdf1_14 , debdf1_15 , debdf1_16 , debdf1_17 , debdf1_18 , debdf1_19 ,debdf1_20 , debdf1_21 , debdf1_22 , debdf1_23 , debdf1_24;
%
%***FIRST EXECUTABLE STATEMENT  INTYD
iflag = 0;
if( k<0 || k>debdf1_20 )
%
iflag = -1;
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;
else;
tp = debdf1_8 - debdf1_7.*(1.0e0+100.0e0.*debdf1_9);
if((t-tp).*(t-debdf1_8)>0.0e0 )
iflag = -2;
else;
%
s =(t-debdf1_8)./debdf1_4;
ic = 1;
if( k~=0 )
jj1 = fix(debdf1_15 - k);
for jj = jj1 : debdf1_20;
ic = fix(ic.*jj);
end; jj = fix(debdf1_20+1);
end;
c = ic;
for i = 1 : debdf1_19;
dky(i) = c.*yh(i,debdf1_15);
end; i = fix(debdf1_19+1);
if( k~=debdf1_20 )
jb2 = fix(debdf1_20 - k);
for jb = 1 : jb2;
j = fix(debdf1_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 : debdf1_19;
dky(i) = c.*yh(i,jp1) + s.*dky(i);
end; i = fix(debdf1_19+1);
end; jb = fix(jb2+1);
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 = debdf1_4.^(-k);
for i = 1 : debdf1_19;
dky(i) = r.*dky(i);
end; i = fix(debdf1_19+1);
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;
%----------------------- end OF SUBROUTINE INTYD -----------------------
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 INVIT