Poincare map of Lorenz equations (homework)

So I am working on a homework assignment and the first part was to plot the lorenz equations on a 3d map. Which I believe I did. Then I needed to modify the code for every time the trajectory crosses the plane z = 27, plot the corresponding point in the xy-plane. And do this for the time interval [30, 70]. Thats where my code gets messed up. This is what I have:

t(1)=0; 
x(1)=1;
y(1)=1;
z(1)=1;
sigma=10;   
rho=28;
beta=(8/3);
h=0.01;   
t=0:h:70;
f=@(t,x,y,z) sigma*(y-x);   
g=@(t,x,y,z) x*rho-x.*z-y;
p=@(t,x,y,z) x.*y-beta*z;
for i=1:(length(t)-1) 
    k1=f(t(i),x(i),y(i),z(i));
    l1=g(t(i),x(i),y(i),z(i));
    m1=p(t(i),x(i),y(i),z(i));
    k2=f(t(i)+h/2,(x(i)+0.5*k1*h),(y(i)+(0.5*l1*h)),(z(i)+(0.5*m1*h)));     
    l2=g(t(i)+h/2,(x(i)+0.5*k1*h),(y(i)+(0.5*l1*h)),(z(i)+(0.5*m1*h)));
    m2=p(t(i)+h/2,(x(i)+0.5*k1*h),(y(i)+(0.5*l1*h)),(z(i)+(0.5*m1*h)));
    k3=f(t(i)+h/2,(x(i)+0.5*k2*h),(y(i)+(0.5*l2*h)),(z(i)+(0.5*m2*h)));
    l3=g(t(i)+h/2,(x(i)+0.5*k2*h),(y(i)+(0.5*l2*h)),(z(i)+(0.5*m2*h)));
    m3=p(t(i)+h/2,(x(i)+0.5*k2*h),(y(i)+(0.5*l2*h)),(z(i)+(0.5*m2*h)));
    k4=f(t(i)+h,(x(i)+k3*h),(y(i)+l3*h),(z(i)+m3*h));
    l4=g(t(i)+h,(x(i)+k3*h),(y(i)+l3*h),(z(i)+m3*h));
    m4=p(t(i)+h,(x(i)+k3*h),(y(i)+l3*h),(z(i)+m3*h));
    x(i+1) = x(i) + h*(k1 +2*k2  +2*k3   +k4)/6; 
    y(i+1) = y(i) + h*(l1  +2*l2   +2*l3    +l4)/6;
    z(i+1) = z(i) + h*(m1+2*m2 +2*m3  +m4)/6;
end
figure(1)
plot3(x,y,z)

for i=30:70 ****************This is where the code gets messed up**************
    for z=27
        x(i)=x;
        y(i)=y;
    end
    figure(2)
    plot(x(i),y(i))
end

My thought process was to use a for loop first with the time interval condition then with the condition that z=27. Then take x(i) which was solved in the runge kutta code and as the value of i gets plugged in when the z value is 27, it would give an x coordinate and y coordinate for when the plane is crossed at z=27. But doing this loop messes up the code to where I cant even plot figure 1 anymore. I feel I am using the for loop wrong, but I thought doing something like this would work. Any hints or documentations to read to help me get on the right track would be greatly appreciated.

Answer 1

[Image Analyst]

Not sure what this means or what your intent is:

x(i)=x;
y(i)=y;

x and y are 7001 long vectors. And you're trying to stuff the whole vector into just one single one of their elements.

If you want to find indexes of x, y, and z where z is 27, then create a mask for where z is in that tight range:

mask = z > 26.9 & z < 27;
% Then extract the x and y values:
x27 = x(mask);
y27 = y(mask);
% Then plot those
plot(x27, y27, 'b-');