program lorenzeuler

      double precision t,x,y,z,dt,tfinal,k(2),l(2),m(2)
      double precision sigma,rho,beta
      integer n,steps

      dt = 0.001
      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 *, f1(x,y,z,sigma,rho,beta)
         k(1) = f1(x,y,z,sigma,rho,beta)*dt
         l(1) = f2(x,y,z,sigma,rho,beta)*dt
         m(1) = f3(x,y,z,sigma,rho,beta)*dt
         k(2) = f1(x+k(1),y+l(1),z+m(1),sigma,rho,beta)*dt
         l(2) = f2(x+k(1),y+l(1),z+m(1),sigma,rho,beta)*dt
         m(2) = f3(x+k(1),y+l(1),z+m(1),sigma,rho,beta)*dt
         t = t + dt
         x = x + 0.5*(k(1) + k(2))
         y = y + 0.5*(l(1) + l(2))
         z = z + 0.5*(m(1) + m(2))

 10   continue

      end


      double precision function f1(x,y,z,sigma,rho,beta)
      double precision x,y,z,sigma,rho,beta
      f1 = sigma*(y-x)
      return
      end

      double precision function f2(x,y,z,sigma,rho,beta)
      double precision x,y,z,sigma,rho,beta
      f2 = x*(rho-z)
      return
      end

      double precision function f3(x,y,z,sigma,rho,beta)
      double precision x,y,z,sigma,rho,beta
      f3 = x*y-beta*z
      return
      end