Event Detection on a circle

I have an ODE trevelling over the surface area.

The ODE starts on the circle however I want the plot to end once it hits the circle again.

Below is the code I am currently using and the photo shows the plot I am currently getting.

The event detection part of the code is the function right at the bottom but I can't seem to get it working.

Any help would be great thanks!

clc

tspan = [0:0.05:80]; %Range for independent var. i.e. time
r = 2;
centreX0 = 0;
centreY0 = 0;
theta = 1.08;
Ptheta = 3;
Pr = -sqrt(2 - 2 * r.^3 .* cos(theta).^3 + 6 .* r.^3 .* cos(theta) .* sin(theta).^2 - (Ptheta.^2 ./ r.^2));
%options = odeset('Events',@MonkeyPlot);

% y0 = [0.03 0.03 0.03 0.03] % Initial Values for x, y, Px, Py
a = r .* cos(theta);
b = r .* sin(theta);
c = Pr .* cos(theta) - (Ptheta ./ r) .* sin(theta);
d = Pr .* sin(theta) + (Ptheta ./ r) .* cos(theta);

y0 = [a b c d]; % Initial Values for x, y, Px, Py

[x, y] = ode45(@ODEdbb,tspan,y0);

[m,n] = meshgrid(-3:.05:3,-3:.05:3);
contour(m,n,m.^3-3*n.^2.*m);

hold on

viscircles([centreX0,centreY0,],r);

axis square

% plot solution of ODE
x_ode = y(:,1);
y_ode = y(:,2);
plot(x_ode,y_ode)
axis([-3 3 -3 3])
hold off

function f = ODEdbb(A, B)
    x=B(1);
    y=B(2);
    Px=B(3);
    Py=B(4);
    
    %Differential Equations
    dPxdt = 3 * y^2 -3 * x^2;
    dPydt = 6 * x * y;
    dxdt = Px;
    dydt = Py;
    f = [dxdt; dydt; dPxdt; dPydt];
end

function[value,isterminal,direction] = MonkeyPlot(A,Y)
   
    value = (Y(1)^2 + Y(2)^2) - r^2;
    isterminal = 1;
    direction = 1;
    
end

Answer 1

[KSSV]

clc

tspan = [0:0.05:80]; %Range for independent var. i.e. time
r = 2;
centreX0 = 0;
centreY0 = 0;
theta = 1.08;
Ptheta = 3;
Pr = -sqrt(2 - 2 * r.^3 .* cos(theta).^3 + 6 .* r.^3 .* cos(theta) .* sin(theta).^2 - (Ptheta.^2 ./ r.^2));
%options = odeset('Events',@MonkeyPlot);

% y0 = [0.03 0.03 0.03 0.03] % Initial Values for x, y, Px, Py
a = r .* cos(theta);
b = r .* sin(theta);
c = Pr .* cos(theta) - (Ptheta ./ r) .* sin(theta);
d = Pr .* sin(theta) + (Ptheta ./ r) .* cos(theta);

y0 = [a b c d]; % Initial Values for x, y, Px, Py

[x, y] = ode45(@ODEdbb,tspan,y0);

Warning: Failure at t=3.877326e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (7.105427e-15) at time t.

% Circle 
viscircles([centreX0,centreY0],r);
hold on

th = linspace(0,2*pi) ;
xc = centreX0+r*cos(th) ;
yc = centreY0+r*sin(th) ; 


[m,n] = meshgrid(-3:.05:3,-3:.05:3);
% Get points insode circle 
idx = inpolygon(m,n,xc,yc) ;
mn = m.^3-3*n.^2.*m ;
mn(~idx) = NaN ; 
contour(m,n,mn);

axis square

% plot solution of ODE
x_ode = y(:,1);
y_ode = y(:,2);
plot(x_ode,y_ode)
axis([-3 3 -3 3])
hold off

function f = ODEdbb(A, B)
    x=B(1);
    y=B(2);
    Px=B(3);
    Py=B(4);
    
    %Differential Equations
    dPxdt = 3 * y^2 -3 * x^2;
    dPydt = 6 * x * y;
    dxdt = Px;
    dydt = Py;
    f = [dxdt; dydt; dPxdt; dPydt];
end