You will need to change your scalar equations to vector equations, and you will need to change your derivative function to a vector function. For the derivative function, use the states x1, x2, x3, x1dot, x2dot, and x3dot in a 6-element vector called y. Then the code would look something like this:
function dy = deriv(t,y,k1,k2,k3,m1,m2,m3)
x1 = y(1);
x2 = y(2);
x3 = y(3);
x1dot = y(4);
x2dot = y(5);
x3dot = y(6);
x1dotdot = ____; % you fill this in
x2dotdot = ____; % you fill this in
x3dotdot = ____; % you fill this in
dy = [x1dot;x2dot;x3dot;x1dotdot;x2dotdot;x3dotdot];
end
You could define your f function in the RK4 routine as (where k1, k2, k3, m1 ,m2, m3 are previously defined):
f = @(t,y) deriv(t,y,k1,k2,k3,m1,m2,m3)
Your RK4 code then needs to be vectorized for the 6-element states. E.g., replace yt(i) with yt(:,i), and replace yt(i-1) with yt(:,i-1). Also the initialization of yt could be this:
yt = zeros(6,ceil(tfinal/h));
yt(:,1) = y0;
