How to plot a Poincare surface of sections

I have a Hamiltonian system with five ODES. I know how to integrtae the equations of motion but I'm not able to plot Poincare surface of sections. I was trying to find it in fileexchange. There is some code here https://it.mathworks.com/matlabcentral/answers/1850203-how-do-i-plot-a-poincare-plot-for-a-non-autonomous-ode-s

But I'm not able to implement it in my code. I'm new learner. Can anyone please help me to plot Poincare surface of sections.

I'm pasting my code here:

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;

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; 
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);

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