How to generalize the RK4 code to solve m number of ODEs in MATLAB?

Question

The following MATLAB code represents the solution of two first order ODEs:


dy1 = @(x,y1,y2) y2; 
dy2 = @(x,y1,y2) 4*x + 10*sin(x) - y1;

%initial values
x0 = pi;
xn = 2*pi;
y1 = 0; % initial value of y1
y2 = 2; % initial value of y2  
h = 0.1; % step size

for x = x0 : h : xn-h
  L1 = h*dy2(x, y1, y2);
  K1 = h*dy1(x, y1, y2);
  L2 = h*dy2(x + h/2, y1 + K1/2, y2 + L1/2);
  K2 = h*dy1(x + h/2, y1 + K1/2, y2 + L1/2);
  L3 = h*dy2(x + h/2, y1 + K2/2, y2 + L2/2);
  K3 = h*dy1(x + h/2, y1 + K2/2, y2 + L2/2);
  L4 = h*dy2(x + h, y1 + K3, y2 + L3);
  K4 = h*dy1(x + h, y1 + K3, y2 + L3);
  y2 = y2 + (L1 + 2*L2 + 2*L3 + L4)/6;
  y1 = y1 + (K1 + 2*K2 + 2*K3 + K4)/6;
  x = x + h;
  fprintf ('%f \t %f\t %f\n',x,y2,y1);
end

How can I make this code general so it could easily manage to solve m number of decomposed first order ODEs (i.e., dy1, dy2, ....., dym)?

Many thanks for help.

Are you looking for something [like this](https://stackoverflow.com/questions/43408704/solve-a-system-of-equations-with-runge-kutta-4-matlab)? — Vicky, Mar 15 '21 at 20:30
@Vicky In fact that, I am looking for something like [https://math.stackexchange.com/questions/721076/help-with-using-the-runge-kutta-4th-order-method-on-a-system-of-2-first-order-od] Specifically, Kai's answer but using MATLAB. — ak ve, Mar 15 '21 at 20:41

James Tursa · Accepted Answer · 2021-03-17T20:58:35.390

Just use state vectors as if you were coding for ode45( ) instead of individual scalar variables for your states. E.g.,

dy = @(x,y) [y(2); 4*x + 10*sin(x) - y(1)];
%initial values
x0 = pi;
xn = 2*pi;
y0 = [0;2]; % initial value of y
h = 0.1; % step size
x = x0 : h : xn-h;
n = numel(x);
y = zeros(numel(y0),n);
y(:,1) = y0;
for i=1:n-1
    K1 = h * dy(x(i)      , y(:,i)       );
    K2 = h * dy(x(i) + h/2, y(:,i) + K1/2);
    K3 = h * dy(x(i) + h/2, y(:,i) + K2/2);
    K4 = h * dy(x(i) + h  , y(:,i) + K3  );
    y(:,i+1) = y(:,i) + (K1 + 2*K2 + 2*K3 + K4)/6;
end

If you ever need to change the number of state variables and differential equations, you simply need to change the function handle and initial conditions. The looping code stays exactly the same. Also, you can pass the function handle, range, and initial conditions to ode45( ) and compare results directly with your results. E.g.,

[x,y] = ode45( dy, [x0 xn], y0 );

How to generalize the RK4 code to solve m number of ODEs in MATLAB?

1 Answers1