How to stop ODE45 using event location

I'm trying to stop ODE45 solver once the solution intersects with the circle of radius 3. Below is the main code with the functions below that. I'm trying to stop the computation of the ODE45 with yV(:,3) = 3(i.e. when it reaches the circle)

import circle.*
import intersections.*
%%%%%%%%%%%%%%%%%%%%%Seting the globla variables%%%%%%%%%%%%%%%%%%%%%%
global X1 Y1 x q R0 G mu H0 cv initial_t initial_dtds initial_r initial_drds initial_theta initial_dthetads initial_phi initial_dphids
global initial_t_CS initial_dtds_CS initial_r_CS initial_drds_CS initial_phi_CS initial_dphids_CS initial_z_CS initial_dzds_CS
s = 0:1:10000;  % time scale in FLRW
q = 0;
%%%%%%%%%%%%%%%%%%%%%Setting the constants %%%%%%%%%%%%%%%%%%%%%%%%%%%%
R0=1;
H0=71000;
cv=3*10^(8);
mu = 1*10^7; 
G = 6.67*10^-11;
%mu = 1; 
%G = 1;

%%%%%%%%%%%%%%%%initial conditions for the RK4 for FLRW%%%%%%%%%%%%%%%%
initial_t    = 1;
initial_dtds = 1;
initial_r    = 5;
initial_drds = -1;
initial_theta  = pi;
initial_dthetads = .05;
initial_phi    = 0;
initial_dphids = 0;
%draw 1st circle at (0,0) radius .5 and get X and Y data
H1=circle([0 0],3,30);
X1=get(H1,'XData');
Y1=get(H1,'YData');
%plotting the circle with dots
plot(X1,Y1, ':','linewidth',1)
hold on
xlim([-1 1]);
ylim([-1 1]);

%%%%%%%%%%%%%%%%the initial condition matrix for FLRW%%%%%%%%%%%%%%%%%%
v0 = [initial_t, initial_dtds, initial_r, initial_drds, ...
      initial_theta, initial_dthetads, initial_phi, initial_dphids];

%%%%%%%%%%%%%%%%%%%%%Solving the rhsflrw equations%%%%%%%%%%%%%%%%%%%%% 
%[s,v] = ode45( @(s,y) rhsflrw(s,y), s, v0);
options = odeset('Events',@events);
[sV, yV] = ode45(@rhsflrw, s, v0,options);
%[x,q]=intersections(X1,Y1,yV(:,3).*cos(yV(:,5)),(yV(:,3).*sin(yV(:,5))),0);

function dxds=rhsflrw(s,y)
global R0 H0 cv initial_t initial_dtds initial_r initial_drds initial_theta initial_dthetads initial_phi initial_dphids

      dxds_1 = y(2);
      dxds_2 = -((2/3)*R0^(2)*((3*H0/(2*cv))^(4/3))*y(1)^(1/3))*(y(4)^(2)+abs(y(3))^(2)*y(6)^(2)+abs(y(3))^(2)*(sin(y(5)))^(2)*y(8)^(2));
      dxds_3 = y(4);
      dxds_4 = -(4/(3*y(1)))*y(2)*y(4)+abs(y(3))*(y(6)^(2)+(sin(y(5)))^(2)*y(8)^(2));
      dxds_5 = y(6);
      dxds_6 = -(2/abs(y(3)))*y(4)*y(6)-(4/(3*y(1)))*y(2)*y(6)+((sin(2*y(5)))*(y(8)^(2))/2);
      dxds_7 = y(8);
      dxds_8 = -2*y(8)*((2/(3*y(1)))*y(2)+(cot(y(5)))*y(6)+(1/abs(y(3)))*y(4));

      dxds=[dxds_1; dxds_2; dxds_3; dxds_4; dxds_5; dxds_6; dxds_7; dxds_8];
  end

function [yV,isterminal,direction] = events(sV,yV);
global X1 Y1 x q R0 G mu H0 cv initial_t initial_dtds initial_r initial_drds initial_theta initial_dthetads initial_phi initial_dphids
yV(:3) = r; % Boundary is r = 3
isterminal = 1;   % Stop the integration
direction = 0;   % You're always approaching from outside
end

Answer 1

[Walter Roberson]

yV(:3) = r is not valid syntax: the ":" operator must either be by itself as an index, or it must have a value on each side of the ":".

You have not defined r anywhere. You have also not defined any other variable with value 3.

Only "global" in the variables you really need, because global is slow.

You should be using something of the form

function [values, isterminal, direction] = events(sV,yV)

where there is one entry in values for each condition. For example you might have

   values = yv(3) - 3

if it just so happened that r is the third variable and you want to stop when it crosses 3.

Further advice: you will find it more reliable to use

values = yv(3) >= 3

as the termination condition. That will become positive when yv(3) reaches 3 from below, and positive is the signal to stop. I have seen people indicating that in practice the ode* routine does not always detect crossings.

[Mr_Turnip_Head]

When I make that edit, I receive an error "Index exceeds matrix dimensions."

[Arnab Joardar]

I tried to follow up your answer on my problem of a bouncing rod. The picture below shows it on the top-left corner. Ground is at Z=-1. The blue line is irrelevant. The red line is the rod.

On the command Window, I am printing the values of the variable 'values' (same variable as in your case). It gives a 0 when the inequality is false; meaning that integration must cease. As you can see, despite 4 occurrences of '0' as one of its values, the 'Event' function did not terminate the integration of the equation of motion. The energy plot on the bottom-left corner shows that the rod had momentarily ceased to move from 0.45 sec to around 5.5 sec. The variable 'Record C' at the bottom shows the instances when collision was recorded between the rod and the ground. I record it after I get a non-empty variable 'IE', which indicates that a collision event has been detected. It is important that I terminate it immediately when the values returned by 'values' is zero.

Can someone please advise me as to how do I resolve this issue of the ODE45 not terminating when explicitly told to? 'isterminal' is 1 and 'direction' is 0.

Thank you in advance.