program rabbit_fox
      implicit double precision (a-h,o-z)
      dimension y(2),aux(2),ak(2,4)
      fun1(yrab,yfox)=1.2d0*yrab-0.6d0*yrab*yfox
      fun2(yrab,yfox)=-0.8d0*yfox+0.3d0*yrab*yfox
c
      time=0.d0
      deltat=0.1
c
      neqs=2
      y(1)=10
      y(2)=2 
      write(*,'(2x,f5.1,10(1x,f10.5))')time,(y(i),i=1,neqs)
c
      do ktimestep=1,500
c-step 1:
         ak(1,1) = deltat*fun1(y(1),y(2))
         ak(2,1) = deltat*fun2(y(1),y(2))
         aux(1) = y(1)+0.5d0*ak(1,1)
         aux(2) = y(2)+0.5d0*ak(2,1)
c-step 2:
         ak(1,2) = deltat*fun1(aux(1),aux(2))
         ak(2,2) = deltat*fun2(aux(1),aux(2))
         aux(1) = y(1)+0.5d0*ak(1,2)
         aux(2) = y(2)+0.5d0*ak(2,2)
c-step 3:
         ak(1,3) = deltat*fun1(aux(1),aux(2))
         ak(2,3) = deltat*fun2(aux(1),aux(2))
         aux(1) = y(1)+ak(1,3)
         aux(2) = y(2)+ak(2,3)
c-step 4:
         ak(1,4) = deltat*fun1(aux(1),aux(2))
         ak(2,4) = deltat*fun2(aux(1),aux(2))
c-synthesis:
      y(1)=y(1)+(ak(1,1)+2.d0*ak(1,2)+2.d0*ak(1,3)+ak(1,4))/6.d0
      y(2)=y(2)+(ak(2,1)+2.d0*ak(2,2)+2.d0*ak(2,3)+ak(2,4))/6.d0
      time=time+deltat
      write(*,'(2x,f5.1,10(1x,f10.5))')time,(y(i),i=1,neqs)
      enddo
c
      end