#### Introduction Numerical methods for solving systems of ordinary differential equations (ODEs) extend scalar techniques to handle multiple coupled equations simultaneously. Such systems are commonly expressed in vector form as $ \frac{d\mathbf{y}}{dt} = \mathbf{f}(t, \mathbf{y}), \quad \mathbf{y}(t_0) = \mathbf{y}_0, $ where $\mathbf{y}(t) \in \mathbb{R}^m$ is the vector of dependent variables, and $\mathbf{f}(t, \mathbf{y})$ is a vector-valued function representing the system of derivatives. The Fourth-Order Runge-Kutta (RK4) method is a widely used, high-accuracy numerical integration technique that evaluates the slope at four points within each time step to achieve fourth-order accuracy. For systems, the RK4 update is given by $ \begin{aligned} \mathbf{k}_1 &= h \mathbf{f}(t_i, \mathbf{y}_i), \\ \mathbf{k}_2 &= h \mathbf{f}\left(t_i + \frac{h}{2}, \mathbf{y}_i + \frac{\mathbf{k}_1}{2}\right), \\ \mathbf{k}_3 &= h \mathbf{f}\left(t_i + \frac{h}{2}, \mathbf{y}_i + \frac{\mathbf{k}_2}{2}\right), \\ \mathbf{k}_4 &= h \mathbf{f}(t_i + h, \mathbf{y}_i + \mathbf{k}_3), \\ \mathbf{y}_{i+1} &= \mathbf{y}_i + \frac{1}{6}(\mathbf{k}_1 + 2\mathbf{k}_2 + 2\mathbf{k}_3 + \mathbf{k}_4), \end{aligned} $ where $h$ is the step size and $t_i = t_0 + i h$. This vectorized RK4 method advances all components of the system simultaneously, combining multiple slope evaluations to provide a highly accurate and stable solution. It balances computational cost and precision, making it a standard choice for numerically solving systems of ODEs. The provided code implements this method for systems, where the function $\mathbf{f}$ returns a row vector of derivatives for each dependent variable. The solution matrix $\mathbf{y}$ stores the approximations for all variables at each time step, with each column representing one component of the system. #### Codes ``` MATLAB function [t, y] = vectorRK4(f, range, iv, h) % Define the time vector from range(1) to range(2) with step size h t = range(1):h:range(2); N = length(t); % Number of time steps % Initialize the solution matrix y with N rows and length(iv) columns y = zeros(N, length(iv)); y(1, :) = iv; % Set the initial condition for i = 1:N-1 % Compute the four slopes k1, k2, k3, and k4 for the RK4 method % k1 is the slope at the beginning of the interval k1 = h * f(t(i), y(i, :)); % k2 is the slope at the midpoint, using k1 for the intermediate value k2 = h * f(t(i) + h/2, y(i, :) + k1/2); % k3 is another slope at the midpoint, using k2 for the intermediate value k3 = h * f(t(i) + h/2, y(i, :) + k2/2); % k4 is the slope at the end of the interval, using k3 for the intermediate value k4 = h * f(t(i) + h, y(i, :) + k3); % Update the solution for the next time step using the weighted average of the slopes y(i+1, :) = y(i, :) + (1/6) * (k1 + 2*k2 + 2*k3 + k4); end end ``` ```python import numpy as np def vector_rk4(f, range_vals, iv, h): """ Vectorized Fourth-Order Runge-Kutta method for systems of ODEs. Parameters: f : function Vector-valued function f(t, y) returning derivatives as a 1D array or list. range_vals : tuple or list (t_start, t_end) integration interval. iv : array-like Initial values for each dependent variable (vector). h : float Step size. Returns: t : numpy.ndarray 1D array of time points. y : numpy.ndarray 2D array of solution values; each column corresponds to one dependent variable. """ t = np.arange(range_vals[0], range_vals[1] + h, h) N = len(t) m = len(iv) y = np.zeros((N, m)) y[0, :] = iv for i in range(N - 1): k1 = h * np.array(f(t[i], y[i, :])) k2 = h * np.array(f(t[i] + h/2, y[i, :] + k1/2)) k3 = h * np.array(f(t[i] + h/2, y[i, :] + k2/2)) k4 = h * np.array(f(t[i] + h, y[i, :] + k3)) y[i + 1, :] = y[i, :] + (k1 + 2*k2 + 2*k3 + k4) / 6 return t, y ```