program rk4sys
      implicit none

      double precision xeuler(2),t,dt,xexact,f,xrk(2),k(4),m(4)
      double precision f1,f2
      open(13,file="data.d")
      dt = 0.125
      xeuler(1) = 1
      xeuler(2) = -1
      xrk(1) = 1
      xrk(2) = -1

      do 10 t=0,10,dt
         xexact = 51.0/47.0*cos(sqrt(1.5)*t) - sqrt(2.0/3.0)*
     c        sin(sqrt(1.5)*t) - 4.0/47.0*cos(5.0*t)
         write(13,11) t,xeuler(1),xexact,xrk(1)
         k(1) = f1(xrk(1), xrk(2),t)*dt
         m(1) = f2(xrk(1), xrk(2),t)*dt

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

         k(3) = f1(xrk(1) + 0.5*k(2), xrk(2) + 0.5*m(2),t + 0.5*dt)*dt
         m(3) = f2(xrk(1) + 0.5*k(2), xrk(2) + 0.5*m(2),t + 0.5*dt)*dt

         k(4) = f1(xrk(1) + k(3), xrk(2) + m(3),t + dt)*dt
         m(4) = f2(xrk(1) + k(3), xrk(2) + m(3),t + dt)*dt

         xrk(1) = xrk(1) + (k(1) + 2.0*k(2) + 2.0*k(3) + k(4))/6.0
         xrk(2) = xrk(2) + (m(1) + 2.0*m(2) + 2.0*m(3) + m(4))/6.0

         xeuler(1) = xeuler(1) + f1(xeuler(1),xeuler(2),t)*dt
         xeuler(2) = xeuler(2) + f2(xeuler(1),xeuler(2),t)*dt
 10   continue
 11   format(4(f15.7))
      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 = 0.5*(4.0*cos(5.0*t)-3.0*x1)
      return
      end