如何在下面的程序中加一个庞加莱截面呢？？谢谢

水风而动

#include <stdio.h>
#include <math.h>
#define N 10000
#define h 0.01
double f(double x,double y)
{
    return(10*(y-x));
}
double p(double x,double y,double z ,double c)
{
  
    return(c*x-x*z-y);
}
double q(double x,double y,double z)
{
    return(x*y-8/3*z);
}
double  runge_kutta(double x0,double y0,double z0)
{  
    double x,y,z,c;
    double d1,d2,d3,d4,k1,k2,k3,k4,f1,f2,f3,f4;
    int i,j;
       x=x0;
    y=y0;
       z=z0;
    FILE *fp;
    fp=fopen("lorenz分岔4.txt","w");
   
    for(j=1;j<=500;j++)
   {

c=1+j;
   
   for(i=1;i<=N;i++)
      {   
       d1=f(x,y);
       k1=p(x,y,z,c);
       f1=q(x,y,z);
      
       d2=f(x+h*d1/2,y+h*k1/2);
       k2=p(x+h*d1/2,y+h*k1/2,z+h*f1/2,c);
       f2=q(x+h*d1/2,y+h*k1/2,z+h*f1/2);
      
       d3=f(x+h*d2/2,y+h*k2/2);
       k3=p(x+h*d2/2,y+h*k2/2,z+h*f2/2,c);
       f3=q(x+h*d2/2,y+h*k2/2,z+h*f2/2);
      
       d4=f(x+h*d3,y+h*k3);
       k4=p(x+h*d3,y+h*k3,z+h*f3,c);
       f4=q(x+h*d3,y+h*k3,z+h*f3);
      
       x=x+h*(d1+2*d2+2*d3+d4)/6;
       y=y+h*(k1+2*k2+2*k3+k4)/6;
       z=z+h*(f1+2*f2+2*f3+f4)/6;
   
        if(i>9500)
        {
      
        fprintf(fp," %8f  , %8f \n",c,x);
        printf("  %8f ,  %8f \n",c,x);
   
        }
     }
  }
}
main()
{  
double x0=-10,y0=0,z0=20;
   
   runge_kutta(x0,y0,z0);  
}

水风而动

这么大的一个论坛居然没有会的吗？