[x,xout,yout,ypout,neq,kold,phi,psi]=ddatrp(x,xout,yout,ypout,neq,kold,phi,psi);

function [x,xout,yout,ypout,neq,kold,phi,psi]=ddatrp(x,xout,yout,ypout,neq,kold,phi,psi);
%***BEGIN PROLOGUE  DDATRP
%***SUBSIDIARY
%***PURPOSE  Interpolation routine for DDASSL.
%***LIBRARY   SLATEC (DASSL)
%***TYPE      doubleprecision (SDATRP-S, DDATRP-D)
%***AUTHOR  Petzold, Linda R., (LLNL)
%***DESCRIPTION
%-----------------------------------------------------------------------
%     THE METHODS IN SUBROUTINE DDASTP USE POLYNOMIALS
%     TO APPROXIMATE THE SOLUTION. DDATRP APPROXIMATES THE
%     SOLUTION AND ITS DERIVATIVE AT TIME XOUT BY EVALUATING
%     ONE OF THESE POLYNOMIALS, AND ITS DERIVATIVE,THERE.
%     INFORMATION DEFINING THIS POLYNOMIAL IS PASSED FROM
%     DDASTP, SO DDATRP CANNOT BE USED ALONE.
%
%     THE PARAMETERS ARE:
%     X     THE CURRENT TIME IN THE INTEGRATION.
%     XOUT  THE TIME AT WHICH THE SOLUTION IS DESIRED
%     YOUT  THE INTERPOLATED APPROXIMATION TO Y AT XOUT
%           (THIS IS OUTPUT)
%     YPOUT THE INTERPOLATED APPROXIMATION TO YPRIME AT XOUT
%           (THIS IS OUTPUT)
%     NEQ   NUMBER OF EQUATIONS
%     KOLD  ORDER USED ON LAST SUCCESSFUL STEP
%     PHI   ARRAY OF SCALED DIVIDED DIFFERENCES OF Y
%     PSI   ARRAY OF PAST STEPSIZE HISTORY
%-----------------------------------------------------------------------
%***ROUTINES CALLED  (NONE)
%***REVISION HISTORY  (YYMMDD)
%   830315  DATE WRITTEN
%   901009  Finished conversion to SLATEC 4.0 format (F.N.Fritsch)
%   901019  Merged changes made by C. Ulrich with SLATEC 4.0 format.
%   901026  Added explicit declarations for all variables and minor
%           cosmetic changes to prologue.  (FNF)
%***end PROLOGUE  DDATRP
%
persistent c d gamma i j koldp1 temp1 ; 

yout_shape=size(yout);yout=reshape(yout,1,[]);
ypout_shape=size(ypout);ypout=reshape(ypout,1,[]);
phi_shape=size(phi);phi=reshape([phi(:).',zeros(1,ceil(numel(phi)./prod([neq])).*prod([neq])-numel(phi))],neq,[]);
psi_shape=size(psi);psi=reshape(psi,1,[]);
%
if isempty(i), i=0; end;
if isempty(j), j=0; end;
if isempty(koldp1), koldp1=0; end;
if isempty(c), c=0; end;
if isempty(d), d=0; end;
if isempty(gamma), gamma=0; end;
if isempty(temp1), temp1=0; end;
%
%***FIRST EXECUTABLE STATEMENT  DDATRP
koldp1 = fix(kold + 1);
temp1 = xout - x;
for i = 1 : neq;
yout(i) = phi(i,1);
ypout(i) = 0.0d0;
end; i = fix(neq+1);
c = 1.0d0;
d = 0.0d0;
gamma = temp1./psi(1);
for j = 2 : koldp1;
d = d.*gamma + c./psi(j-1);
c = c.*gamma;
gamma =(temp1+psi(j-1))./psi(j);
for i = 1 : neq;
yout(i) = yout(i) + c.*phi(i,j);
ypout(i) = ypout(i) + d.*phi(i,j);
end; i = fix(neq+1);
end; j = fix(koldp1+1);
%
%------end OF SUBROUTINE DDATRP------
yout_shape=zeros(yout_shape);yout_shape(:)=yout(1:numel(yout_shape));yout=yout_shape;
ypout_shape=zeros(ypout_shape);ypout_shape(:)=ypout(1:numel(ypout_shape));ypout=ypout_shape;
phi_shape=zeros(phi_shape);phi_shape(:)=phi(1:numel(phi_shape));phi=phi_shape;
psi_shape=zeros(psi_shape);psi_shape(:)=psi(1:numel(psi_shape));psi=psi_shape;
end
%DECK DDAWS