Large system of ordinary differential equations I am trying to find a large system (>20) of coupled ordinary differential equations in order to approximate them numerically on the  computer and check the efficiency and effectiveness of various steppers. Can someone point me to such a system? Thank you in advance
 A: If you don't care about the type of system (i.e. whether it's stiff or not) you could consider the [Heat equation][1].
$$\frac{\partial u }{\partial t} -\lambda \Delta^2 u = f$$
with $ \Delta u =\nabla^2 u = \sum_i \frac{\partial^2 u}{\partial x_i^2}$  and $\lambda > 0$.
Let us consider a 2D domain of length 1. We define $\Omega =[0,1]\times [0,1]$. If we want to approximate our domain with $n$ elements in each direction, we get $n^2$ elements in total. Therefore our equidistant stepsize  $h=\frac{1}{n+1}$ if we exclude (count only inner points) the boundary. Let us assume $u=0$ on $\partial \Omega$, so we have homogenous boundary condition, and $u(t=0) = 0$ as our initial condition We approximate $\Delta u$ by $h^2\Delta u(i\cdot h, j\cdot h) \approx u_{i-1,j}+u_{i+1,j}-4u_{ij} +u_{i,j+1}+u_{i,j-1}$. With that we get:
\begin{align}
&u(t=0) = 0 \quad \text{as the initial condition}\\
&u(\partial \Omega) = 0 \quad \text{as our boundary condition}\\
&\frac{\partial}{\partial t} u_{i,j} - \frac{\lambda}{h^2} (u_{i-1,j}+u_{i+1,j}-4u_{ij} +u_{i,j+1}+u_{i,j-1}) = f_{ij}\\
\Leftrightarrow & \frac{\partial}{\partial t} u_{i,j} = \frac{\lambda}{h^2} (u_{i-1,j}+u_{i+1,j}-4u_{ij} +u_{i,j+1}+u_{i,j-1}) + f_{ij}\\
\end{align}
Now we have to set up our grid. Let in the following $G$ denote the inner grid points, that means that $G$ does not include the boundary. We will enumerate our elements as follows:
\begin{align}
G = \begin{pmatrix}
1 & 2 & 3 & ... & n \\
n+1 & n+2 & n+3 &... & 2n\\
... & ...& ... & ... &... &\\
... & ...& ... & ... &n^2
\end{pmatrix}
\end{align}
This should present the numbering of $u$, so $u(2\cdot h, 3 \cdot h) = u_{n+3}$. Of course the same holds for $f$.
If we rewrite our ODE as $\tilde u_t =\frac{\lambda}{h^2} A \tilde u+F$  we get 
\begin{align}
\tilde u &= (u_1,u_2,...,u_{n^2})^T\\
F &= (f_1,f_2,...,f_{n^2})^T \\
A &= \begin{pmatrix}
\tilde{A} & I & 0 & 0 & ... & 0\\
I & \tilde{A} & I & 0 & ... & 0\\
0 &I & \tilde{A} & I &  ... & 0 \\
... & ...& ...& ...& ...& ...\\
... & ...& ...& ...& \tilde{A} & I 
\end{pmatrix} \in \mathbb{R}^{n^2 \times n^2}
\end{align}
with 
\begin{align}
\tilde{A} &= \text{tridiag}(1,-4,1) \in \mathbb{R}^{n \times n}\\
 I &= \text{diag(1)} \in \mathbb{R}^{n \times n}\\
 0  &\in \mathbb{R}^{n \times n}
\end{align}
This allows you to set up large systems easily, which has a physicall meaningfull background. For a simple example, you could choose $F$ to be zero everywhere except on one side of $\Omega$, heat should diffuse through the domain.
Note, that your system tends to get large very soon, as $A \in \mathbb{R}^{n^2\times n^2}$. Also note, that this is a sparse matrix, so there might be a method which has an advantage if it can handle sparse matrices.
