Program ODE1

      Real*8 dx, x0, y0, z0, t0, xi, yi, zi, ti, xf, yf, zf, tf, tmax
      Integer*4 i
      External dx, dy, dz
      open (unit=7,file="ode1d.dat",status='unknown')
c* initial conditions
      x0 = 0.0
      y0 = 2.0
      z0 = 0.0
      t0 = 0.00
c* "time" step and max time
      dt =   0.01
      tmax = 100.00

      xi = x0
      yi = y0
      zi = z0
      ti = t0
      do while(ti.le.tmax)
         write (7,*) xi, yi,zi
         tf = ti + dt
         call  rk4_1(dx,ti,xi,yi,zi,tf,xf)
         xi = xf
         call  rk4_2(dy,ti,xi,yi,zi,tf,yf)
         yi = yf
         call  rk4_3(dz,ti,xi,yi,zi,tf,zf)
         zi = zf
	   ti = tf
      end do
c100   format(2e12.4)
      stop
      end

c-------------------------------------------
      Function dx(t,x,y,z)
      Real*8 dx, x, y, z, t, a
      a=10
      dx = -a*(x-y)
      return
      end
c-------------------------------------------
      Function dy(t,x,y,z)
      Real*8 dy,x, y, z, t, r
      r=28
      dy = -x*z+r*x-y
      return
      end

c-------------------------------------------
      Function dz(t,x,y,z)
      Real*8 dz,x, y, z, t, b
      b=8/3
      dz = x*y-b*z
      return
      end
c-------------------------------------------
      Subroutine rk4_1(dx,ti,xi,yi,zi,tf,xf)

      Real*8 dx,ti,xi,yi,zi,tf,xf
      Real*8 h,k1,k2,k3,k4

      h  = tf-ti

      k1 = h*dx(ti,xi,yi,zi)
      k2 = h*dx(ti+h/2.0,xi+k1/2.0,yi,zi)
      k3 = h*dx(ti+h/2.0,xi+k2/2.0,yi,zi)
      k4 = h*dx(ti+h,xi+k3,yi,zi)

      xf = xi + (k1 + 2.0*(k2+k3) + k4)/6.0
      Return
      End
c-------------------------------------------
      Subroutine rk4_2(dy,ti,xi,yi,zi,tf,yf)

      Real*8 dy,ti,xi,yi,zi,tf,yf
      Real*8 h,k1,k2,k3,k4

      h  = tf-ti

      k1 = h*dy(ti,xi,yi,zi)
      k2 = h*dy(ti+h/2.0,xi,yi+k1/2.0,zi)
      k3 = h*dy(ti+h/2.0,xi,yi+k2/2.0,zi)
      k4 = h*dy(ti+h,xi,yi+k3,zi)

      yf = yi + (k1 + 2.0*(k2+k3) + k4)/6.0
      Return
      End
c-------------------------------------------
      Subroutine rk4_3(dz,ti,xi,yi,zi,tf,zf)

      Real*8 dz,ti,xi,yi,zi,tf,zf
      Real*8 h,k1,k2,k3,k4

      h  = tf-ti

      k1 = h*dz(ti,xi,yi,zi)
      k2 = h*dz(ti+h/2.0,xi,yi,zi+k1/2.0)
      k3 = h*dz(ti+h/2.0,xi,yi,zi+k2/2.0)
      k4 = h*dz(ti+h,xi,yi,zi+k3)

      zf = zi + (k1 + 2.0*(k2+k3) + k4)/6.0
      Return
      End