Problem with coding for Poincare surface of sections

Im trying to plot Poincare surface of sections. I know how to integrate equations of motion, but having difficulty to define the conditions and collecting points to plot Poincare sections. I have plotted something, but it's wrong - I have checked my plot with Mathematica plot.

I'm pasting my code here. Could anyone please help?? I shall be highly thankful for you.

clear All

clc

EE=1.38; LL=40; BB = 1;

%here r = y(1), theta = y(2), \phi = y(3), p_r = y(4), p_\theta = y(5)

f =@(y) [y(4)*(1-2/y(1)); y(5)*(1/(y(1))^2); -BB + (LL*(csc(y(2)))^2)/(y(1)^2);

((1/(y(1))^3)*((y(5))^2))-(1/(y(1))^2)*((y(4))^2)+(1/2)*((-2*(EE^2))/(((-1+(2/y(1)))^2)*(y(1))^2)+(2*((LL - BB*(y(1)^2)*(sin(y(2))^2))^2)*(csc(y(2)))^2)/((y(1))^3) + (4*BB*(LL - BB*(y(1)^2)*(sin(y(2))^2)))/(y(1)));

(1/2)*(4*BB*(cot(y(2)))*(LL - BB*(y(1)^2)*(sin(y(2))^2)) + (2/(y(1))^2)*((LL - BB*(y(1)^2)*(sin(y(2))^2))^2)*(cot(y(2)))*(csc(y(2)))*(csc(y(2))))];

h = 0.1; N = 1000-1; y(:,1) = [8; 1.06; 0; 0; 0]; t(1)=0.0;

y0=y(:,1);

% Hamiltonian at initial conditions

He=(1/2)+(1/2)*(1-(2)/(y0(1)))*((y0(4))^2)+((y0(5))^2)/(2*(y0(1))^2)+(1/2)*((EE^2)/(-1+(2)/(y0(1)))+((((LL - BB*(y0(1)^2)*(sin(y0(2))^2))^2)*(csc(y0(2)))*(csc(y0(2))))/(y0(1)^2)));

% Butcher Table

a21=1/2; a31=0; a32=1/2; a41=0; a42=0; a43=1; a44=0;t=0;

c1=0; c2=1/2; c3=1/2; c4=1;

b1=1/6; b2=1/3; b3=1/3; b4=1/6;

for n=1:N

k1=y(:,n);

k2=y(:,n)+h*a21*(f(k1));

k3=y(:,n)+h*a31*(f(k1))+h*a32*(f(k2));

k4=y(:,n)+h*a41*(f(k1))+h*a42*(f(k2))+h*a43*(f(k3));

y(:,n+1)=y(:,n)+h*b1*(f(k1))+h*b2*(f(k2))+h*b3*(f(k3))+h*b4*(f(k4));

t(n+1)=t(n)+h;

end

% Hamiltonian at output value

Hn=(1./2)+(1./2).*(1-2./y(1,:)).*(y(4,:).^2)+(y(5,:).^2)./(2.*(y(1,:).^2))+(1./2).*((EE^2)./(-1+(2)./(y(1,:)))+(((LL - BB.*(y(1,:).^2).*(sin(y(2,:)).^2)).^2).*(csc(y(2,:))).*(csc(y(2,:))))./((y(1,:)).^2));

error= abs(He-Hn);

y;

% Extract relevant data for Poincaré surface

r = y(1, :); theta = y(2, :); phi = y(3, :); p_r = y(4, :); p_theta = y(5, :);

% Define Poincaré section condition

threshold =pi/2; % Threshold value for phi

% Find indices where the condition is satisfied

poincareIndices = find(diff(theta) > 0 & theta(1:end-1) < threshold);

% Parameters for Poincaré section

sectionTheta = pi/2; % Poincaré section angle (choose the desired value)

sectionTolerance = 1e-6; % Tolerance for sectionTheta condition

% Extract data corresponding to Poincaré section

poincareR = r(poincareIndices);

poincareTheta = theta(poincareIndices);

poincarephi = phi(poincareIndices);

poincarep_r = p_r(poincareIndices);

poincarep_theta = p_theta(poincareIndices);

% Plot Poincaré surface

figure;

scatter(poincareR,poincarep_r,'.', 'MarkerEdgeColor', 'g');

xlabel('R');

ylabel('p_r');

title('Plot b/w r and p_r');

%plot3(cos (y(3,:)) .*y (1,:) .*sin(y(2,:)),sin (y(3,:)) .*y (1,:) .*sin(y(2,:)),y (1,:) .*cos(y(2,:)))