Why are 2 events being listed instead of only one ?

Hi I am solving a second order differential equation dealing with expansion and collapse of a gas bubble in water. Events is used to detect the end of collapse and expansion (In eventsFuncn1A).

The end of collapse is detected using the condition that -dR/dt = 0 and with slope of velocity being negative, i.e dR/dt is positive. (event 1) The end of expansion is detected using the condition that dR/dt = 0 and with slope of velocity being negative. (event 2)

On running the program, 2 events are being listed. Shouldn't only 1 event be listed because the stop condition in 'eventsFuncn1A' is that the program should stop at the end of either collapse or expansion.

The output also shows that the 2 events listed are the same (Time 6.260857e-005, Radius 4.707202e-004, Velcoity -4.913847e-013)

I can't understand why length(tevent) should give 2 events. Since I need the no. of events for further processing, i am unable to proceed.

Will someone please explain this.

TIA

*The programs are given below :*
% MAIN SCRIPT
clear all;
close all;
clc;

tspan = [0 1e-3]; z0 = [0.1e-3 ; 0]  % bubl size < resonance
T0=293 ;

options = odeset('stats', 'on', 'outputfcn', 'odeplot', 'Events', ...
@eventsFuncn1A);

gamma = 1.33;
sig = 0.0725; % Surface tension

pvp=0  % For AIR bubble

patm = 1*10^5;

pa = 2.4*10^5;  

w = 10000*2*pi;  % frequency = 10 kHz

mu = 1e-3; % Viscosity (kg/ms) 

rho = 1000;
R0 = z0(1) ; 
p(1) = w; p(2) = pa; p(3) = patm; p(4)= gamma; p(5) = sig;
p(6) = pvp; p(7) = mu; p(8) = R0; p(9) = rho;

times = [];
solns = [];

[t z tevent zevent indexevent] = ode23s(@zdotA, tspan, z0, options, p); 

format short e;
nsteps = length(t);
times = [times; t(1:nsteps)]
solns = [solns; z(1:nsteps,:)]

nevents = length(tevent);
fprintf('\nnumber of events  %d', nevents)

for n=1:nevents

     fprintf('\n EVENT NUMBER  %d', n)
     fprintf('\n INDEX EVENT %d', indexevent(n))
     fprintf('\n Time %d', tevent(n)), fprintf('\n Radius %d', zevent(n, 

1)),
   fprintf('\n Velcoity %d', zevent(n,2))

if (indexevent(n)==1 )
display('CHANGING GAMMA - now starting expansion')
    gamma = 1;
p(4) = gamma

elseif (indexevent(n)==2)
      display('CHANGING GAMMA - now starting compression')  
    gamma = 1.33;
    p(4) = gamma;
  end
end

z0 = z(nsteps, :);  % changes initial condition for next expansion / 

contraction
z0(2) = -z0(2);   
tspan(1) = t(nsteps);

x = z(:, 1); y = z(:, 2);
T = T0*(x/R0).^(3*(gamma-1));
% Pressure at interface
pg0 = p(3) + 2*p(5)/p(8) - p(6);  
%pg0 = 0;   %(NOLTINGK & NEPPIRAS)- Fig. 2-CASE OF VOID (no gas inside) 

pg = pg0 * (R0./z(:,1)).^(3*gamma);

pr = p(6) + pg - (2*sig./x) - 4*mu*y./x;
whpr = rho*y*1450;

subplot(4,2,1); plot(t,x)
grid on
%axis([0 1e-3 1.7e-3  2.5e-3]) % LEIGHTON Fig 4.6 only
xlabel('time'); ylabel('Radius');

subplot(4,2,2); plot(t,y)
xlabel('time'); ylabel('Velocity');
disp('Max jet velocity ')
max(abs(y))

%subplot(4,2,3); plot(y.*t/R0, x./R0);
%xlabel('Ut/R0'); ylabel('R/R0');
subplot(4,2,3); plot(t, x./R0);
xlabel('t'); ylabel('R/R0');
%subplot(4,2,4); plot(y.*t/0.01, x./0.01)
%xlabel('Ut/R0'); ylabel('R/R0');
%subplot(4,2,4); plot(w*t, x./R0)
%xlabel('w*t'); ylabel('R/R0');

subplot(4,2,4); plot(t,pr,'--')
grid on
xlabel('time'); ylabel('pressure in liquid');
disp('Max pressure in liquid')
max(pr)

subplot(4,2,5); plot(t,p(6)+pg)
grid on
xlabel('time'); ylabel('pressure inside bubble');
disp('Max pressure inside bubble ')
max(p(6)+pg)

%subplot(4,2,5); plot(t,whpr)
%xlabel('time'); ylabel('water hammer pressure');

subplot(4,2,6); plot(t,T)
xlabel('time'); ylabel('Internal Temp, K');

disp('Max water hammer pressure')
max(whpr)

*%%%%%%%%%%%
% Events Function*
function [eventvalue, stopthecalc eventdirection] = eventsFuncn1A(t, z, p)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
R = z(1);
Rdot = z(2);
eventvalue = [-Rdot Rdot];
stopthecalc = [1 1];
eventdirection = [1 -1];

end

%%%%%%%%%%%%%%
% Function returning derivative
function [rdot] = zdotA(t, z, p)
% function [rdot, B] = zdot4(t, z, p)  HOW TO PASS B OUTSIDE?
w = p(1); pa = p(2); p0 = p(3); gamma = p(4); sig = p(5); 
pvp = p(6); mu = p(7); R0 = p(8); rho = p(9);

A = p0 - pa*sin(w*t) ;

pg0 = p0 - pvp + 2*sig/R0; % pg0 = patm - pvp + 2*sig/R

B = pg0 * (R0 / z(1))^(3*gamma);

C = 4*mu*z(2)/z(1);
D = 2*sig/z(1); % ST of water from EXERC IN ADVANC COMP MECH

rr2dot = ((pvp+B-A-C-D)/rho - 3/2*z(2)^2)/z(1);

rdot = [z(2); rr2dot];

end