[h,k,n,nq,t,tout,yh,y]=ddntp(h,k,n,nq,t,tout,yh,y);

function [h,k,n,nq,t,tout,yh,y]=ddntp(h,k,n,nq,t,tout,yh,y);
%***BEGIN PROLOGUE  DDNTP
%***SUBSIDIARY
%***PURPOSE  Subroutine DDNTP interpolates the K-th derivative of Y at
%            TOUT, using the data in the YH array.  If K has a value
%            greater than NQ, the NQ-th derivative is calculated.
%***LIBRARY   SLATEC (SDRIVE)
%***TYPE      doubleprecision (SDNTP-S, DDNTP-D, CDNTP-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
%***ROUTINES CALLED  (NONE)
%***REVISION HISTORY  (YYMMDD)
%   790601  DATE WRITTEN
%   900329  Initial submission to SLATEC.
%***end PROLOGUE  DDNTP
persistent factor i j jj kk kused r ; 

if isempty(i), i=0; end;
if isempty(j), j=0; end;
if isempty(jj), jj=0; end;
if isempty(kk), kk=0; end;
if isempty(kused), kused=0; end;
if isempty(factor), factor=0; end;
if isempty(r), r=0; end;
y_shape=size(y);y=reshape(y,1,[]);
yh_shape=size(yh);yh=reshape([yh(:).',zeros(1,ceil(numel(yh)./prod([n])).*prod([n])-numel(yh))],n,[]);
%***FIRST EXECUTABLE STATEMENT  DDNTP
if( k==0 )
for i = 1 : n;
y(i) = yh(i,nq+1);
end; i = fix(n+1);
r =((tout-t)./h);
for jj = 1 : nq;
j = fix(nq + 1 - jj);
for i = 1 : n;
y(i) = yh(i,j) + r.*y(i);
end; i = fix(n+1);
end; jj = fix(nq+1);
else;
kused = fix(min(k,nq));
factor = 1.0d0;
for kk = 1 : kused;
factor = factor.*(nq+1-kk);
end; kk = fix(kused+1);
for i = 1 : n;
y(i) = factor.*yh(i,nq+1);
end; i = fix(n+1);
r =((tout-t)./h);
for jj = kused + 1 : nq;
j = fix(kused + 1 + nq - jj);
factor = 1.0d0;
for kk = 1 : kused;
factor = factor.*(j-kk);
end; kk = fix(kused+1);
for i = 1 : n;
y(i) = factor.*yh(i,j) + r.*y(i);
end; i = fix(n+1);
end; jj = fix(nq+1);
for i = 1 : n;
y(i) = y(i).*h.^(-kused);
end; i = fix(n+1);
end;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
yh_shape=zeros(yh_shape);yh_shape(:)=yh(1:numel(yh_shape));yh=yh_shape;
end
%DECK DDOGLG