Need help to stop the ode when v(2) is zero using ode events. I have tried in so many ways but couldn't get result.

%%Data for desnity with respect to depth 
z = [2 3 5 7 10 15 20 25 30 40 50 60 70 80 90 100 125 150 160 175 200 225 250 275 300 325 350 375 400];
rho = [17.2731684 17.1649375 21.43455647 22.4140625 23.86332207 24.3746967 24.70487685 24.6003125 24.8933125 25.42772826 26.03220776 26.439625 26.8151875 26.86830797 27.1949375 27.34406944 27.5551875 27.728625 27.23423729 27.88542857 27.752249049 28.1025 28.2415 28.37 28.05366667 28.6565 28.7755 28.898 29.013]; 
v0=[0.8;0.001;0.1;]; % Initial values 
% Creating a matrix 
zsol = [];
v1sol = [];
v2sol = [];
v3sol = [];
for i=1:numel(rho) - 1
rho0=17.1;
g=9.8;
zspan = [z(i) z(i+1)];
Nsquare = (g/rho0)*(rho(i+1)-rho(i))/(z(i+1)-z(i));
options=odeset('RelTol',1e-4,'AbsTol',1e-6,'Events',@events);
[t,v] = ode45(@(t,v)rhs(t,v,Nsquare), zspan, v0,options);
v0 = [v(end,1) ; v(end,2) ; v(end,3)];
zsol =  [zsol;t];
v1sol = [v1sol;v(:,1)];
v2sol = [v2sol;v(:,2)];
v3sol = [v3sol;v(:,3)];
end
plot(v1sol,zsol,'r',v2sol,zsol,'g',v3sol,zsol,'b')
xlabel('Width,b and vertical velocity,w')
ylabel('Height, z')
grid on ;
function parameters=rhs(~,v,Nsquare)
        alpha=0.116;
            db= 4*alpha*v(2).^2-v(1).*v(3)./2*v(2).^2;
            dw= v(1).*v(3)-4*alpha*v(2).^4+v(1).*v(2).^2.*v(3)./2*v(1).*v(2).^3;
            dgmark= (-2*Nsquare*v(1).*v(2)^4-v(1).*v(3)^2+4*alpha*v(2).^4.*v(3)-v(1).*v(2).^2.*v(3).^2-8*alpha*v(2).^3.*v(3)+2*v(3).^2.*v(1).*v(2))./2*v(1).*v(2).^4; 
            parameters=[db;dw;dgmark];
end

[Walter Roberson]

You did not post code for events and the MATLAB function named events would not be relevant to this matter ?

[Dereje]

yes your are right. I posted the original code. here is the event codes I used in the program.

 options=odeset('RelTol',1e-4,'AbsTol',1e-6,'Events',@events);
 [t,v] = ode45(@(t,v)rhs(t,v,Nsquare), zspan, v0,options);
  %
%
function [value, isterminal, direction] = events(t,v)
   value = v(2);
   isterminal = 1;
   direction = 0;
end

Answer 1

[Jan]

But your works sufficiently already. Insert this test:

[t,v] = ode45(@(t,v)rhs(t,v,Nsquare), zspan, v0,options);
if t(end) ~= zspan
  disp('Event function has stopped the integration');
end

You will see, that the event function works as expected. At i=14 the 0 is reached by v(2). So what does "couldn't get result" mean?

Note that you have some other troubles:

 Warning: Failure at t=1.802503e+02.  Unable to meet integration tolerances
 without reducing the step size below the smallest value allowed
 (4.547474e-13) at time t.

It seems like there is a pole. Then the integration cannot work reliably.

[Dereje]

By adjusting the initial values of vo possible to get a good result. But still I am not able to see v(2) zero at i = 14(like you said) in the result graph. How could I able to see. Also no matter what intital values I used for vo, It can not stop.

[Jan]

@Dereje: I cannot follow you. The event function finds a zero crossing for i=14. I cannot imagine, why you are not "able to see it". Where do you try to check this?

i = 14;   % Instead of the FOR loop
rho0=17.1;
g=9.8;
zspan = [z(i) z(i+1)];
Nsquare = (g/rho0)*(rho(i+1)-rho(i))/(z(i+1)-z(i));
options=odeset('RelTol',1e-4,'AbsTol',1e-6,'Events',@events);
[t,v] = ode45(@(t,v)rhs(t,v,Nsquare), zspan, v0,options);
plot(t, v);

 Also no matter what intital values I used for vo, It can not stop.

Please be clear: what cannot "stop" and what do you mean by "stopping"? The integration stops, if v(2) is 0 in your code already, so what do you want instead?

[Dereje]

Hi Jan, Thanks a lot! Now I understand. I was mixing how the integration stops and the choice of initial values(The graph of my original codes, v(1),v(2) and v(3) varies greatly with small change in initial values.)

[Dereje]

Also When I say It can not stop what I meant was: let's say v(2) is zero at z=80. Though in the result graph there is values for z greater than 80 (In principle I would have to see results only for z 0 to 80).

[Dereje]

From the graph you can see that v(2) is zero around 80 and you can also see that for z larger than 80.