Efficient solution of ODE with variable dimensions

I am solving a system of ODEs, that consists of only one equation at t = 0. As the solution progresses in time, more and more equations get added to the system at specific events. Eventually, the system can possibly have tens of thousands of dimensions. The max. number of dimensions are unknown a priori.

The way I am solving the problem is to create a very large system of equations, most of which are unused most of the time. This is very inefficient.

Is there a way such a problem could be solved more efficiently - such that dimensions are added and removed as required?

The jacobian of the system cannot be evaluated in closed form.

Any help will be appreciated!