[x,y,xout,yout,ypout,neqn,kold,phi,ivc,iv,kgi,gi,alpha,og,ow,ox,oy]=sintrp(x,y,xout,yout,ypout,neqn,kold,phi,ivc,iv,kgi,gi,alpha,og,ow,ox,oy);

function [x,y,xout,yout,ypout,neqn,kold,phi,ivc,iv,kgi,gi,alpha,og,ow,ox,oy]=sintrp(x,y,xout,yout,ypout,neqn,kold,phi,ivc,iv,kgi,gi,alpha,og,ow,ox,oy);
persistent alp c g gamma gdi gdif h hi hmu i iq iw j jq kp1 kp2 l m rmu sigma temp1 temp2 temp3 w xi xim1 xiq ; 

if isempty(alp), alp=0; end;
if isempty(c), c=zeros(1,13); end;
if isempty(g), g=zeros(1,13); end;
if isempty(gamma), gamma=0; end;
if isempty(gdi), gdi=0; end;
if isempty(gdif), gdif=0; end;
if isempty(h), h=0; end;
if isempty(hi), hi=0; end;
if isempty(hmu), hmu=0; end;
if isempty(rmu), rmu=0; end;
if isempty(sigma), sigma=0; end;
if isempty(temp1), temp1=0; end;
if isempty(temp2), temp2=0; end;
if isempty(temp3), temp3=0; end;
if isempty(w), w=zeros(1,13); end;
if isempty(xi), xi=0; end;
if isempty(xim1), xim1=0; end;
if isempty(xiq), xiq=0; end;
if isempty(i), i=0; end;
if isempty(iq), iq=0; end;
if isempty(iw), iw=0; end;
if isempty(j), j=0; end;
if isempty(jq), jq=0; end;
if isempty(kp1), kp1=0; end;
if isempty(kp2), kp2=0; end;
if isempty(l), l=0; end;
if isempty(m), m=0; end;
%***BEGIN PROLOGUE  SINTRP
%***PURPOSE  Approximate the solution at XOUT by evaluating the
%            polynomial computed in STEPS at XOUT.  Must be used in
%            conjunction with STEPS.
%***LIBRARY   SLATEC (DEPAC)
%***CATEGORY  I1A1B
%***TYPE      SINGLE PRECISION (SINTRP-S, DINTP-D)
%***KEYWORDS  ADAMS METHOD, DEPAC, INITIAL VALUE PROBLEMS, ODE,
%             ORDINARY DIFFERENTIAL EQUATIONS, PREDICTOR-CORRECTOR,
%             SMOOTH INTERPOLANT
%***AUTHOR  Watts, H. A., (SNLA)
%***DESCRIPTION
%
%   The methods in subroutine  STEPS  approximate the solution near  X
%   by a polynomial.  Subroutine  SINTRP  approximates the solution at
%   XOUT  by evaluating the polynomial there.  Information defining this
%   polynomial is passed from  STEPS  so  SINTRP  cannot be used alone.
%
%   subroutine STEPS is completely explained and documented in the text,
%   'Computer Solution of Ordinary Differential Equations, the Initial
%   Value Problem'  by L. F. Shampine and M. K. Gordon.
%
%   Input to SINTRP --
%
%   The user provides storage in the calling program for the arrays in
%   the call list
%      DIMENSION Y(NEQN),YOUT(NEQN),YPOUT(NEQN),PHI(NEQN,16),OY(NEQN)
%                AND ALPHA(12),OG(13),OW(12),GI(11),IV(10)
%   and defines
%      XOUT -- point at which solution is desired.
%   The remaining parameters are defined in  STEPS  and passed to
%   SINTRP  from that subroutine
%
%   Output from  SINTRP --
%
%      YOUT(*) -- solution at  XOUT
%      YPOUT(*) -- derivative of solution at  XOUT
%   The remaining parameters are returned unaltered from their input
%   values.  Integration with  STEPS  may be continued.
%
%***REFERENCES  H. A. Watts, A smoother interpolant for DE/STEP, INTRP
%                 II, Report SAND84-0293, Sandia Laboratories, 1984.
%***ROUTINES CALLED  (NONE)
%***REVISION HISTORY  (YYMMDD)
%   840201  DATE WRITTEN
%   890831  Modified array declarations.  (WRB)
%   890831  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   920501  Reformatted the REFERENCES section.  (WRB)
%***end PROLOGUE  SINTRP
%
y_shape=size(y);y=reshape(y,1,[]);
yout_shape=size(yout);yout=reshape(yout,1,[]);
ypout_shape=size(ypout);ypout=reshape(ypout,1,[]);
phi_orig=phi;phi_shape=[neqn,16];phi=reshape([phi_orig(1:min(prod(phi_shape),numel(phi_orig))),zeros(1,max(0,prod(phi_shape)-numel(phi_orig)))],phi_shape);
oy_shape=size(oy);oy=reshape(oy,1,[]);
%
%***FIRST EXECUTABLE STATEMENT  SINTRP
kp1 = fix(kold + 1);
kp2 = fix(kold + 2);
%
hi = xout - ox;
h = x - ox;
xi = hi./h;
xim1 = xi - 1.;
%
%   INITIALIZE W(*) FOR COMPUTING G(*)
%
xiq = xi;
for iq = 1 : kp1;
xiq = xi.*xiq;
temp1 = iq.*(iq+1);
w(iq) = xiq./temp1;
end; iq = fix(kp1+1);
%
%   COMPUTE THE DOUBLE INTEGRAL TERM GDI
%
if( kold<=kgi )
gdi = gi(kold);
else;
if( ivc>0 )
iw = fix(iv(ivc));
gdi = ow(iw);
m = fix(kold - iw + 3);
else;
gdi = 1.0./temp1;
m = 2;
end;
if( m<=kold )
for i = m : kold;
gdi = ow(kp2-i) - alpha(i).*gdi;
end; i = fix(kold+1);
end;
end;
%
%   COMPUTE G(*) AND C(*)
%
g(1) = xi;
g(2) = 0.5.*xi.*xi;
c(1) = 1.0;
c(2) = xi;
if( kold>=2 )
for i = 2 : kold;
alp = alpha(i);
gamma = 1.0 + xim1.*alp;
l = fix(kp2 - i);
for jq = 1 : l;
w(jq) = gamma.*w(jq) - alp.*w(jq+1);
end; jq = fix(l+1);
g(i+1) = w(1);
c(i+1) = gamma.*c(i);
end; i = fix(kold+1);
end;
%
%   DEFINE INTERPOLATION PARAMETERS
%
sigma =(w(2)-xim1.*w(1))./gdi;
rmu = xim1.*c(kp1)./gdi;
hmu = rmu./h;
%
%   INTERPOLATE FOR THE SOLUTION -- YOUT
%   AND FOR THE DERIVATIVE OF THE SOLUTION -- YPOUT
%
for l = 1 : neqn;
yout(l) = 0.0;
ypout(l) = 0.0;
end; l = fix(neqn+1);
for j = 1 : kold;
i = fix(kp2 - j);
gdif = og(i) - og(i-1);
temp2 =(g(i)-g(i-1)) - sigma.*gdif;
temp3 =(c(i)-c(i-1)) + rmu.*gdif;
for l = 1 : neqn;
yout(l) = yout(l) + temp2.*phi(l,i);
ypout(l) = ypout(l) + temp3.*phi(l,i);
end; l = fix(neqn+1);
end; j = fix(kold+1);
for l = 1 : neqn;
yout(l) =((1.0-sigma).*oy(l)+sigma.*y(l))+ h.*(yout(l)+(g(1)-sigma.*og(1)).*phi(l,1));
ypout(l) = hmu.*(oy(l)-y(l))+(ypout(l)+(c(1)+rmu.*og(1)).*phi(l,1));
end; l = fix(neqn+1);
%
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
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_orig(1:prod(phi_shape))=phi;phi=phi_orig;
oy_shape=zeros(oy_shape);oy_shape(:)=oy(1:numel(oy_shape));oy=oy_shape;
end
%DECK SIR