program lorenz

c     program to determine an approximate numerical solution to the
c     Lorenz equations.

      double precision x,y,z,sigma,beta,rho,t,dt
      common /coef/ sigma,beta,rho
      double precision k(4),l(4),m(4)

      dt=0.05

      x = 0.1
      y = 0.1
      z = 0.1 
      
      sigma = 10.0
      beta = 8.0/3.0
      rho = 28.0

      open(unit=44,file="lorenz.d")

      do 10 t=0,50,dt

         k(1) = f(x,y,z,t)*dt
         l(1) = g(x,y,z,t)*dt
         m(1) = h(x,y,z,t)*dt

         k(2) = f(x+k(1)/2.0,y+l(1)/2.0,z+m(1)/2.0,t+dt/2.0)*dt
         l(2) = g(x+k(1)/2.0,y+l(1)/2.0,z+m(1)/2.0,t+dt/2.0)*dt
         m(2) = h(x+k(1)/2.0,y+l(1)/2.0,z+m(1)/2.0,t+dt/2.0)*dt

         k(3) = f(x+k(2)/2.0,y+l(2)/2.0,z+m(2)/2.0,t+dt/2.0)*dt
         l(3) = g(x+k(2)/2.0,y+l(2)/2.0,z+m(2)/2.0,t+dt/2.0)*dt
         m(3) = h(x+k(2)/2.0,y+l(2)/2.0,z+m(2)/2.0,t+dt/2.0)*dt

         k(4) = f(x+k(3), y+l(3), z+m(3), t+dt)*dt
         l(4) = g(x+k(3), y+l(3), z+m(3), t+dt)*dt
         m(4) = h(x+k(3), y+l(3), z+m(3), t+dt)*dt

         x = x + (k(1) + 2.0*k(2) + 2.0*k(3) + k(4))/6.0
         y = y + (l(1) + 2.0*l(2) + 2.0*l(3) + l(4))/6.0
         z = z + (m(1) + 2.0*m(2) + 2.0*m(3) + m(4))/6.0

         write(44,*) t,x,y,z

 10   continue 

      stop
      end

      double precision function f(x,y,z,t)
      double precision x,y,z,t,sigma,beta,rho
      common /coef/ sigma,beta,rho
      f = sigma*(y - x)
      return
      end

      double precision function g(x,y,z,t)
      double precision x,y,z,t,sigma,beta,rho
      common /coef/ sigma,beta,rho
      g = x*(rho - z)-y
      return
      end

      double precision function h(x,y,z,t)
      double precision x,y,z,t,sigma,beta,rho
      common /coef/ sigma,beta,rho
      h = x*y - beta*z
      return
      end