program duffing4rk

      double precision x(2),t,dt,tfinish,k(4),l(4)
      integer n,steps

      dt = 0.01
      tfinish = 50
      steps = tfinish/dt+1
      print *, steps
      t = 0
      x(1) = 1
      x(2) = 1

      open(unit=13,file="data.d")

      do 10 n=1,steps
         print *, t
         write(13,*) t,x(1),x(2)
         t = t + dt      

         k(1) = f1(x(1),          x(2),          t)*dt
         l(1) = f2(x(1),          x(2),          t)*dt

         k(2) = f1(x(1)+0.5*k(1), x(2)+0.5*l(1), t+0.5*dt)*dt
         l(2) = f2(x(1)+0.5*k(1), x(2)+0.5*l(1), t+0.5*dt)*dt

         k(3) = f1(x(1)+0.5*k(2), x(2)+0.5*l(2), t+0.5*dt)*dt
         l(3) = f2(x(1)+0.5*k(2), x(2)+0.5*l(2), t+0.5*dt)*dt

         k(4) = f1(x(1)+k(3),     x(2)+l(3),     t+dt)*dt
         l(4) = f2(x(1)+k(3),     x(2)+l(3),     t+dt)*dt

         x(1) = x(1) + (k(1)+2*k(2)+2*k(3)+k(4))/6.0D0
         x(2) = x(2) + (l(1)+2*l(2)+2*l(3)+l(4))/6.0D0
 10   continue
      end

      double precision function f1(x1,x2,t)
      double precision x1,x2,t
      f1 = x2
      return
      end

      double precision function f2(x1,x2,t)
      double precision x1,x2,t
      f2 =  cos(t) - 0.1*x2 + x1 - x1**3
      return
      end