[f,neq,t,y,h,yp,f1,f2,f3,f4,f5,ys,rpar,ipar]=defehl(f,neq,t,y,h,yp,f1,f2,f3,f4,f5,ys,rpar,ipar);

function [f,neq,t,y,h,yp,f1,f2,f3,f4,f5,ys,rpar,ipar]=defehl(f,neq,t,y,h,yp,f1,f2,f3,f4,f5,ys,rpar,ipar);
persistent ch k ; 

if isempty(ch), ch=0; end;
if isempty(k), k=0; end;
%***BEGIN PROLOGUE  DEFEHL
%***SUBSIDIARY
%***PURPOSE  Subsidiary to DERKF
%***LIBRARY   SLATEC
%***TYPE      SINGLE PRECISION (DEFEHL-S, DFEHL-D)
%***AUTHOR  Watts, H. A., (SNLA)
%***DESCRIPTION
%
%     Fehlberg Fourth-Fifth order Runge-Kutta Method
% **********************************************************************
%
%    DEFEHL integrates a system of NEQ first order
%    ordinary differential equations of the form
%               dU/DX = F(X,U)
%    over one step when the vector Y(*) of initial values for U(*) and
%    the vector YP(*) of initial derivatives, satisfying  YP = F(T,Y),
%    are given at the starting point X=T.
%
%    DEFEHL advances the solution over the fixed step H and returns
%    the fifth order (sixth order accurate locally) solution
%    approximation at T+H in the array YS(*).
%    F1,---,F5 are arrays of dimension NEQ which are needed
%    for internal storage.
%    The formulas have been grouped to control loss of significance.
%    DEFEHL should be called with an H not smaller than 13 units of
%    roundoff in T so that the various independent arguments can be
%    distinguished.
%
%    This subroutine has been written with all variables and statement
%    numbers entirely compatible with DERKFS. For greater efficiency,
%    the call to DEFEHL can be replaced by the module beginning with
%    line 222 and extending to the last line just before the return
%    statement.
%
% **********************************************************************
%
%***SEE ALSO  DERKF
%***ROUTINES CALLED  (NONE)
%***REVISION HISTORY  (YYMMDD)
%   800501  DATE WRITTEN
%   890831  Modified array declarations.  (WRB)
%   891009  Removed unreferenced statement label.  (WRB)
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900328  Added TYPE section.  (WRB)
%   910722  Updated AUTHOR section.  (ALS)
%***end PROLOGUE  DEFEHL
%
%
y_shape=size(y);y=reshape(y,1,[]);
yp_shape=size(yp);yp=reshape(yp,1,[]);
f1_shape=size(f1);f1=reshape(f1,1,[]);
f2_shape=size(f2);f2=reshape(f2,1,[]);
f3_shape=size(f3);f3=reshape(f3,1,[]);
f4_shape=size(f4);f4=reshape(f4,1,[]);
f5_shape=size(f5);f5=reshape(f5,1,[]);
ys_shape=size(ys);ys=reshape(ys,1,[]);
rpar_shape=size(rpar);rpar=reshape(rpar,1,[]);
ipar_shape=size(ipar);ipar=reshape(ipar,1,[]);
%
%***FIRST EXECUTABLE STATEMENT  DEFEHL
ch = h./4.;
for k = 1 : neq;
ys(k) = y(k) + ch.*yp(k);
end; k = fix(neq+1);
[dumvar1,ys,f1,rpar,ipar]=f(t+ch,ys,f1,rpar,ipar);
%
ch = 3..*h./32.;
for k = 1 : neq;
ys(k) = y(k) + ch.*(yp(k)+3..*f1(k));
end; k = fix(neq+1);
[dumvar1,ys,f2,rpar,ipar]=f(t+3..*h./8.,ys,f2,rpar,ipar);
%
ch = h./2197.;
for k = 1 : neq;
ys(k) = y(k) + ch.*(1932..*yp(k)+(7296..*f2(k)-7200..*f1(k)));
end; k = fix(neq+1);
[dumvar1,ys,f3,rpar,ipar]=f(t+12..*h./13.,ys,f3,rpar,ipar);
%
ch = h./4104.;
for k = 1 : neq;
ys(k) = y(k)+ ch.*((8341..*yp(k)-845..*f3(k))+(29440..*f2(k)-32832..*f1(k)));
end; k = fix(neq+1);
[dumvar1,ys,f4,rpar,ipar]=f(t+h,ys,f4,rpar,ipar);
%
ch = h./20520.;
for k = 1 : neq;
ys(k) = y(k) + ch.*((-6080..*yp(k)+(9295..*f3(k)-5643..*f4(k)))+(41040..*f1(k)-28352..*f2(k)));
end; k = fix(neq+1);
[dumvar1,ys,f5,rpar,ipar]=f(t+h./2.,ys,f5,rpar,ipar);
%
%     COMPUTE APPROXIMATE SOLUTION AT T+H
%
ch = h./7618050.;
for k = 1 : neq;
ys(k) = y(k)+ ch.*((902880..*yp(k)+(3855735..*f3(k)-1371249..*f4(k)))+(3953664..*f2(k)+277020..*f5(k)));
end; k = fix(neq+1);
%
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
yp_shape=zeros(yp_shape);yp_shape(:)=yp(1:numel(yp_shape));yp=yp_shape;
f1_shape=zeros(f1_shape);f1_shape(:)=f1(1:numel(f1_shape));f1=f1_shape;
f2_shape=zeros(f2_shape);f2_shape(:)=f2(1:numel(f2_shape));f2=f2_shape;
f3_shape=zeros(f3_shape);f3_shape(:)=f3(1:numel(f3_shape));f3=f3_shape;
f4_shape=zeros(f4_shape);f4_shape(:)=f4(1:numel(f4_shape));f4=f4_shape;
f5_shape=zeros(f5_shape);f5_shape(:)=f5(1:numel(f5_shape));f5=f5_shape;
ys_shape=zeros(ys_shape);ys_shape(:)=ys(1:numel(ys_shape));ys=ys_shape;
rpar_shape=zeros(rpar_shape);rpar_shape(:)=rpar(1:numel(rpar_shape));rpar=rpar_shape;
ipar_shape=zeros(ipar_shape);ipar_shape(:)=ipar(1:numel(ipar_shape));ipar=ipar_shape;
end
%DECK DEFER