program rkex4

      implicit none

      double precision x1,x2,f1,f2,k(4),m(4),dt,t
      open(13,file="data.d")
      dt = 0.5
      x1 = 1
      x2 = 0

      do 10 t=0,20,dt
         write(13,11) t,x1,x2,x1-cos(t)

         k(1) = f1(x1           ,x2,        t)*dt
         m(1) = f2(x1           ,x2         ,t)*dt

         k(2) = f1(x1 + 0.5*k(1),x2+0.5*m(1),t+0.5*dt)*dt
         m(2) = f2(x1 + 0.5*k(1),x2+0.5*m(1),t+0.5*dt)*dt

         k(3) = f1(x1 + 0.5*k(2),x2+0.5*m(2),t+0.5*dt)*dt
         m(3) = f2(x1 + 0.5*k(2),x2+0.5*m(2),t+0.5*dt)*dt

         k(4) = f1(x1 + k(3)    ,x2 + m(3) ,t + dt)*dt
         m(4) = f2(x1 + k(3),    x2 + m(3),  t + dt)*dt

         x1 = x1 + 1.0/6.0*(k(1) + 2.0*k(2) + 2.0*k(3) + k(4)) 
         x2 = x2 + 1.0/6.0*(m(1) + 2.0*m(2) + 2.0*m(3) + m(4)) 

 10   continue
 11   format(4(f15.6))

      stop
      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 = -x1

      return
      end