real t,x,y,z,dt,tfinal
      real sigma,rho,beta
      real f1,f2,f3
      integer n,steps

      dt = 0.005
      tfinal = 40
      steps = tfinal/dt+1

      sigma = 10
      rho = 10
      beta = 8.0/3.0
      x = 1
      y = 1
      z = 1

      open(unit=13,file="data.d")

      do 10 n=1,steps
         write(13,*) t,x,y,z
         print *, t
         t = t + dt
         f1 = sigma*(y - x)
         f2 = x*(rho - z)-y
         f3 = x*y - beta*z
         x = x + f1*dt + dt**2/2.0*(-sigma*f1 + sigma*f2)
         y = y + f2*dt + dt**2/2.0*((rho-z)*f1 - f2 - x*f3)
         z = z + f3*dt + dt**2/2.0*(y*f1 + x*f2 -beta*f3)

 10   continue

      end