1) Define your state as a column vector, the number of elements being the number of equations that you have (since they are all 1st order equations). If there were higher order equations then we would need more states.
2) Create a function file that calculates the derivative of your state vector at a given time. Each element of this vector will correspond to one of your derivative equations. Use the same signature for this function that the MATLAB ode functions use:
% File my_derivative.m
function dy = my_derivative(t,y)
% put code here to calculate dy as a function of t and y
return
end
E.g., for these equations
y1' = -y2 - y3;
y2' = y1 + (k1*y2);
y3' = k2 + y3*(y1 - k3);
with k1=k2= 0.2 ; k3=5.7
you could have something like this:
% File my_derivative.m
function dy = my_derivative(t,y)
k1 = 0.2;
k2 = 0.2;
k3 = 5.7;
dy = zeros(size(y));
dy(1) = -y(2) - y(3);
dy(2) = y(1) + k1*y(2);
dy(3) = k2 + y(3)*(y(1) - k3);
return
end
Alternatively, if you wanted to pass in the k's as parameters, you could do this:
% File my_derivative.m
function dy = my_derivative(t,y,k1,k2,k3)
dy = zeros(size(y));
dy(1) = -y(2) - y(3);
dy(2) = y(1) + k1*y(2);
dy(3) = k2 + y(3)*(y(1) - k3);
return
end
3) Then in your loop you only have one set of RK4 equations that operates on your state vector, instead of multiple RK4 equations that each operate only on a scalar. Put the iterative results in columns of the result. E.g., a simple outline would be this:
y(:,1) = your initial column state vector
t(1) = the initial time
for i=1:N-1
% Rewrite your derivative equations using (:,i) for indexing instead of (i)
% You will only have y(:,i) here, not y1(:,i) and y2(:,i) etc.
% Be sure to use t(i) in your equations instead of just t (if there is a time dependency)
% Note that you will have only one set of RK4 equations here
y(:,i+1) = y(:,i) + (1/6)*(s1+2*s2+2*s3+s4)*h; % store iterative results in columns of y
t(i+1) = t(i) + h;
end
