      program lorenz
      INTEGER NSTEP,NVAR
      PARAMETER(NVAR=3,NSTEP=6000)
      INTEGER i,j
      REAL xx(NSTEP),x1,x2,y(NVAR,NSTEP),vstart(NVAR)
      REAL sigma,r,b
      COMMON /parameters/ sigma,r,b
      COMMON /path/ xx,y
      EXTERNAL derivs
      x1=1.0
 

      sigma=10.0
      b=8.0/3.0
      r=28.0

      x1=0.0
      x2=60.0
      vstart(1)=0.0
      vstart(2)=0.10
      vstart(3)=0.0

      call rkdumb(vstart,NVAR,x1,x2,NSTEP,derivs)

      open (unit=10,file="lorenz-xy.dat",status="unknown")
      open (unit=11,file="lorenz-xz.dat",status="unknown")
      open (unit=12,file="lorenz-yz.dat",status="unknown")
      open (unit=13,file="lorenz-xyz.dat",status="unknown")

      do 1 i=1,nstep
      write(10,*) y(1,i),y(2,i)
      write(11,*) y(1,i),y(3,i)
      write(12,*) y(2,i),y(3,i)
      write(13,*) y(1,i),y(2,i),y(3,i)
1     continue

      close(10)
      close(11)
      close(12)
      close(13)
      stop
      end

      subroutine derivs(x,y,dydx)

      real x,y(*),dydx(*),sigma,r,b
      common /parameters/ sigma,r,b 

      dydx(1)=-sigma*y(1)+sigma*y(2)
      dydx(2)=-y(1)*y(3)+r*y(1)-y(2)
      dydx(3)=y(1)*y(2)-b*y(3)

      return
      end


      SUBROUTINE rkdumb(vstart,nvar,x1,x2,nstep,derivs)
      INTEGER nstep,nvar,NMAX,NSTPMX
      PARAMETER (NMAX=3,NSTPMX=6000)
      REAL x1,x2,vstart(nvar),xx(NSTPMX),y(NMAX,NSTPMX)
      EXTERNAL derivs
      COMMON /path/ xx,y
CU    USES rk4
      INTEGER i,k
      REAL h,x,dv(NMAX),v(NMAX)
      do 11 i=1,nvar
        v(i)=vstart(i)
        y(i,1)=v(i)
11    continue
      xx(1)=x1
      x=x1
      h=(x2-x1)/nstep
      do 13 k=1,nstep
        call derivs(x,v,dv)
        call rk4(v,dv,nvar,x,h,v,derivs)
        if(x+h.eq.x)pause 'stepsize not significant in rkdumb'
        x=x+h
        xx(k+1)=x
        do 12 i=1,nvar
          y(i,k+1)=v(i)
12      continue
13    continue
      return
      END
