ode23 Integral Limits

I have been attempting to set multiple limits on my system of ODE's. I wrote a script that works, but is slow, and probably not the correct way to do this (attached with working script). This is a simple example of what I am trying to do. In the future I hope to run ten's of thousands of iterations so I am looking for something more efficient.

I attemped to use Event functions to do this by recreating the bouncing ball example. However in that example the ball is only at it's limit for one instance in time. In my case the limit can hold over many time steps. This means I can hit the limit (let's say zero), but the next time step will also stay at the limit (and so on), until the motion is reveresed.

I am using ode23 if the following way.

1) Making the time range only over three time steps (short amount of time). I needed to use 3 time steps because ode23 is a second order method.

%Time Range

dt = 0.05;

t0 = 0.00;

tend = t0 + 2*dt;

2) Setting the Initial Condtions y0.

3) Defining a system of equations for a forced mass on a spring.

function dydt = MassSpring(t,y)

%System of Differential Equations

dydt = [y(2); -0.05*y(1) + 0.1*sin(0.2*pi*t)];

end

4) Running ode23 over the one increment of time (three steps for 2nd order solver).

%Call ODE solver (Only Solve One Time Step)

[t,y] = ode23s(@MassSpring,[t0:dt:tend],y0);

5) Defining upper and lower limits to the travel of the forced mass on a spring.

%Upper and Lower Limits of Integration

a1 = 0.00;

a2 = 0.25;

6) Checking if the solution to one time step of the ODE is between a1 and a2. If yes creating new intial conditions based on ode23 solution.

if and (y(n,1) > a1, y(n-1,1) > a1, y(n,1) < a2, y(n-1,1) < a2)

y0 = [y(n,1); y(n,2)];

7) Checking if the solution is less than the lower limit a1. If it is I change the output of my ode 'y' to lower limits, and update intial condition with limits.

elseif or (y(n,1) < a1, y(n-1,1) < a1)

y0 = [a1; 0.0];

y = [y(n-2,1), a1, a1; y(n-2,2), 0.0, 0.0].';

8) Checking if the solution is greter than the upper limit a2. If it is I change the output of my ode 'y' to upper limits, and updating initial condition with limits.

elseif or (y(n-1,1) > a2, y(n,1) > a2)

y0 = [a2; 0.0];

y = [y(n-2,1), a2, a2; y(n-2,2), 0.0, 0.0].';

9) Creating a list of t and y outputs sucsessivly. This way I can take each ode23 output, run it through my conditinal limit and create a new tout, yout to be plotted.

tout = [tout; t(2:n)];

yout = [yout; y(2:n,:)];

10) I loop on this to my desired solution time. This means I am unloading and loading my ode23 each loop and checking the end of the solution to make sure it stays inbetween it's limits.

[Matthew Hollister]

Can anyone tell me how to this with Event functions?