I’m not sure if you can do exactly what you want, but it is possible to do quite a lot with events. First, this sounds like some sort of numerical calculation of first passage time (a.k.a first hitting time). If these “particles” are stochastic, stop and don’t use ode45 but rather a a method appropriate for SDEs.
As far as I know there is no limit on how many event functions you can have – or rather the dimension of the event function (similar to the dimension of your ODE function) – and their number is not tied to how many ODE equations you have. The events function receives both the current time and the current state vector. You can any or all of the elements of those to create each event. You are correct that more events functions and more complex events will slow down integration. The performance also depends on how often events are detected. If each of your particles reaches the “roof,” as you call it, and triggers just a single event then that won’t be too bad.
In terms of implementation, here a simple example based on Matlab’s ballode example, that simulates N ballistic particles in the vertical alone. There are N non-terminating events to catch the time and speed of each particle as it passes through y = 0. One additional terminating event is added to check if all of the particles have passed through y = 0 (if we knew which particle this would be, as we do here, we could just make that event a terminating one).
function eventsdemo
% Initial conditions for n balls
n = 10;
y0(2*n,1) = 0;
y0(n+1:end) = linspace(20,40,n);
% Specify events function
options = odeset('Events',@(t,y)efun(t,y,n));
% Integrate
[t,y,te,ye,ie] = ode45(@(t,y)f(t,y,n),[0 10],y0,options);
figure;
plot(t,y(:,1:n),'b',te(1:n),ye(sub2ind(size(ye),ie(1:n),(1:n).')),'r.');
function dydt = f(t,y,n)
% Differential equations for ballistic motion
dydt = [y(n+1:end);zeros(n,1)-9.8];
function [value,isterminal,direction] = efun(t,y,n)
% Last event checks that all balls have hit ground and terminates integration
yn = y(1:n);
value = [yn;all(yn < 0)];
zn = zeros(n,1);
isterminal = [zn;1];
direction = [zn-1;1];
In some ways this is slightly inefficient because we keep simulating all N systems even after some of them have passed through y = 0. However, it is simple and the output arrays are rectangular rather than ragged.
I’m not clear what you mean precisely by “coupled together” and needing the particles to “stop.” If you need to do more than just record the the event data, such as change system parameters or change the the differential equation(s) in some other way, then you’ll need to terminate after each event and restart the integration. Look at the ballode example (type edit ballode in the Matlab command window) to see some suggestions on to make this a bit more efficient.