I am looking for a way to share the output of my event function with the main function used by the ode solver.

I have a Monte Carlo simulation which I am trying to make more efficient. For each simulation, random initial conditions are generated and a couple of ode's are solved. Using an event function, the ode's finish computing when certain events occur. Currently, this is done using a for loop and takes 10 hours.

If I look at a simplified example, it seems that it is about 60 times faster to solve all the ode's in one go. That is, I feed the function a matrix of initial conditions (each row representing one set of initial conditions) and get back a 3-dimensional array representing the solution.

But my simplified example doesn't use an event function. It just ends at some arbitrary time. I have been trying (and failing) to figure out how to get my events function to "talk" to my main ode function so that I can turn off the calculation for simulations that are done. Or, more generally, I want to figure out at each time step which of my ode's is still "live" and pass that on to the next time step.

Any help much appreciated. For what it is worth, my code for the simplified example is below.

nsim=1000;
y0=rand(nsim,2);
FuncHandle=(@(t,y)ODEMatrixTester(t,y,nsim));
[t,y]=ode113(FuncHandle,[0 20],y0);

z=zeros(size(t,1),nsim,2);    %Answer comes back as matrix.  This turns it
z(:,:,1)=y(:,1:nsim);         %into a 3-dimensional array.  The front layer
z(:,:,2)=y(:,nsim+1:2*nsim);  %is the first variable, the back the second.
                              %The rows are time steps.  Each column is one
                              %simulation.

function dn = ODEMatrixTester( t, n, nsim )

%This function gives back dn.  Basically, given a vector of values for "n",
%this calculate the rhs of the differential equations. This is fed back 
%to ODE45 (or similar ODE solver) which then knows dn/dt at that point.  It
%moves forward one time step and estimates the changes to n.  It then 
%feeds the new n back to this function.  
%This function is the van der pol equation.  

n=reshape(n,[nsim,2]);  %initial conditions passed as a vector.  
                        %This turns them back into a matrix

dn = [n(:,2),(1-n(:,1).^2).*n(:,2)-n(:,1)];
dn=dn(:);

end

Answer 1

[Josh Meyer]

If you specify a terminal event in the event function, then the integration will halt when the event occurs. You can then use the output returned by the event function (te,ye,ie) to alter the initial conditions and restart the integration.

More information: ODE Event Location

[Cambrian45]

Thank you Joshua. That is very helpful.