1function [t x rec iend] = rk4fixed(t, x, rec, simcons, h, Nstart, xsample)
 2% Matlab implementation of a fixed-step RK4 algorithm.
 3% Written by R.G. Longoria 950129.
 4% This version was revised July 1998 to be compatible with
 5% the Matlab ODE suite.
 6% Usage: >> [T,Y] = rk4fixed('Fcn',[T0 Tf],X0,N);
 7% Inputs to program:
 8% X0 = initial value(s) of the dependent variable.
 9%      Example: X0=[1;1];
10% Tspan = vector, [T0 Tf]
11% T0 = starting value of the independent variable (e.g. time)
12% Tf = ending value of the independent variable
13% h = time step
14% Fcn = the name of the function that defines the ODEs, f(t,x).
15%     Note: User needs to define this function in an m-file:
16%           function yprime=Fcn(T,X)
17% Computed in rk4fixed:
18% h = fixed step size = (Tspan(2)-Tspan(1))/N
19% Output from rk4fixed:
20% Two arrays: X and T over the time period specified.
21%
22% Here is an example of the use of this program from\\
23% the command line:
24%      [T,X]=rk4fixed(`myfunction',[0 0.5],[0.1;0.0],500);
25%      This integrates 2 eqs specified in `myfunction' from
26%      time 0 to 0.5 with steps of h=(0.5-0)/500.  The initial
27%      conditions are given by [0.1;0.0].
28%
29
30%#codegen
31
32iend = 0;
33halfh = 0.5*h;
34N = length(t);
35
36% i_end = Nstart;
37
38td = t(Nstart);
39xd = zeros(size(xsample));
40xd(:) = x(:,Nstart);
41
42for i = (Nstart+1):N
43    [RK1 rec(:, i) events] = equations_of_motion(td, xd, t, x, rec, i, h, simcons);
44    iend = i;
45    if events
46        break
47    end
48    
49    Thalf = td + halfh;
50    Xtemp = xd + halfh*RK1;
51    RK2 = equations_of_motion(Thalf, Xtemp, t, x, rec, i, h, simcons);
52    
53    Xtemp = xd + halfh*RK2;
54    RK3 = equations_of_motion(Thalf, Xtemp, t, x, rec, i, h, simcons);
55    
56    Tfull = td + h;
57    Xtemp = xd + h*RK3;
58    RK4 = equations_of_motion(Tfull, Xtemp, t, x, rec, i, h, simcons);
59    
60    x(:,i) = xd + h*(RK1+2.0*(RK2+RK3)+RK4)/6;
61    
62    xd(:) = x(:,i);
63    td = t(i);  
64end
65