How to solve a set of ODE with discontinuities?

i've been stuck on this problem for a long time now, and i wasn't able to find a solution to this.

the only solution i know of is to use the event function for the solver of the ODE, order it to terminate the integration at the event points, then restart the integration. i've tried implementing this in a while loop, where in each iteration of the loop i itegrate from the previous event point to the next event point. the loop looks something like this:

     ...
    x=[];
    y=[];

while max(x_i)<x_final 
    x_span_i=[x_event_prev,x_final] 
    [x_i,y_i,x_event_i,y_event_i,i_event_i]=ode15s(@D_eqs,x_span_i,y_event_prev,options);
    x=[x;x_i]; 
    y=[y;y_i]; 
    x_event_prev=x_event_i; 
    y_event_prev=y_event_i;
end

where D_eqs is my dy/dx function, and my options and event function are:

tol=10^-9;

      function [value,indicator_terminate,direction]=my_event(x,y) 
      value=[f(y(4))-1,f(y(5))-1,abs(y(4))-tol,abs(y(5))-tol];
      indicator_terminate=[1,1,1,1];
      direction=[0,0,0,0]; 
      end

options=odeset('Events',@my_event,'Refine',20,'RelTol',10^-5,'AbsTol',10^-8);

f is a function defined earlier in the code. to be more precise f is given as follows:

Gamma_thr=10^6;

      function f=f(y) 
      f=(Gamma_thr^-2+y.^2).^(1/2); 
      end

the discotinuities in the equations (i have a set of 5 ODE's) are at:

f(y(4))=1 , f(y(5))=1

the other 2 events (abs(y(4))=tol and abs(y(5))=tol) are just to make the solver stop whenever y(4) or y(5) reaches zero.

my attempt in obtaining a solution doesn't work when i try picking different values for some of my parameters. when this happens i noticed that if i add the line:

x_event_i

inside the while loop, so that it will print its value on each iteration of loop, i get almost right from the start that x_event_i is a vector with 2 values. how can this be that there is more than 1 event point in the iteration, if the itegration is suppost to terminate as soon as one of the events happens?

in addition to that, i get that the solver is pretty much stuck inside the while loop as it prints event points at very very small distance from each other (maybe as much as if every single point is an event point) which obviously makes no sense because y(x) is suppost to evovle with x and therefore the event function f(y(4)) or f(y(5)) should also evovle with x.

i get the feeling that i am missing something here, and i need your help in finding the solution for this numeric/matlab programing problem.

thanks.