1function [t x rec iend] = rk4fixed(t, x, rec, simcons, h, Nstart, xsample)
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32iend = 0;
33halfh = 0.5*h;
34N = length(t);
35
36
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