#### Introduction
Numerical methods for solving systems of ordinary differential equations (ODEs) extend scalar ODE 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 [[Explicit Euler Method]] generalizes naturally to systems by updating the solution vector at each time step using the slope evaluated at the current point:
$
\mathbf{y}_{i+1} = \mathbf{y}_i + h \, \mathbf{f}(t_i, \mathbf{y}_i),
$
where $h$ is the step size, and $t_i = t_0 + i h$ for $i = 0, 1, 2, \ldots$.
This vectorized Euler method approximates the solution by stepping forward in time for all components simultaneously. While straightforward and computationally efficient, it remains only first-order accurate and can suffer from stability issues, especially for stiff systems.
The provided code implements this vectorized Euler method, where the function $\mathbf{f}$ returns a row vector of derivatives corresponding to each dependent variable. The solution matrix $\mathbf{w}$ stores the approximations for all variables at each time step, with each column representing one component of the system.
This approach forms the basis for extending higher-order methods, such as Runge-Kutta schemes, to systems of ODEs, enabling more accurate and stable numerical integration of complex dynamical models.
#### Codes
```MATLAB
function [t,w] = eulerSys(f,range,iv,h)
% input f is a vector valued function containing the system of differential equations
% f takes t and w as inputs, where w is an array containing the dependent variables
% the output f should be organized as a row vector
% input range is a vector of length 2 containing an initial time and a final time
% input iv is a vector containing the initial conditions for each dep variable
% input h is a scalar step size
% output t is a column vector containing the times at which we have solution approximations
% output w is a matrix whose columns are the approximations of the different dep variables
%set up an evenly spaced array for t
t = (range(1):h:range(2))';
% Preallocate the solution matrix w. Use the length of t to determine the number of time steps
% and the length of iv to determine the number of functions/ODEs.
% Make sure the matrix is tall, not wide. That is, each column of w should represent one function.
n=length(t);
m=length(iv);
%input the initial conditions into the first row of w
w = zeros(n,m);
w(1,:)=iv;
%run through Euler's method to fill in the remaining rows of w
for i=1:n-1
w(i+1,:) = w(i,:) + h*f(t(i),w(i,:));
end
```
```python
import numpy as np
def euler_sys(f, range_vals, iv, h):
"""
Explicit Euler 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.
w : 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)
w = np.zeros((n, m))
w[0, :] = iv
for i in range(n - 1):
w[i + 1, :] = w[i, :] + h * np.array(f(t[i], w[i, :]))
return t, w
```