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

-
You can create arbitrarily large linear ODE systems of the form $y'=Ay+F(t)$ for an appropriate sized $A$. Alternatively, you can create non linear ODE systems by choosing the functions yourself. In this manner, you would have a benchmark to test the numerical error of your solver as you know the exact solution(s). –  Daryl Feb 9 '13 at 23:23

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.