program ex

      double precision x(2),t,dt,tfinish,k(4),m(4)

      x(1) = 1
      x(2) = 0

      tfinish=20
      dt= 0.25

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

      do 10 t=0,tfinish,dt
         write(13,*) t,x(1),cos(t),cos(t)-x(1)

         k(1) = f1(x(1),x(2),t)*dt
         k(2) = f1(x(1)+k(1)/2,x(2)+m(1)/2,t+dt/2)*dt
         k(3) = f1(x(1)+k(2)/2,x(2)+m(2)/2,t+dt/2)*dt
         k(4) = f1(x(1)+k(3),x(2)+m(3),t+dt)*dt

         m(1) = f2(x(1),x(2),t)*dt
         m(2) = f2(x(1)+k(1)/2,x(2)+m(1)/2,t+dt/2)*dt
         m(3) = f2(x(1)+k(2)/2,x(2)+m(2)/2,t+dt/2)*dt
         m(4) = f2(x(1)+k(3),x(2)+m(3),t+dt)*dt

         x(1) = x(1) + (k(1) + 2*k(2) + 2*k(3) + k(4))/6.0
         x(2) = x(2) + (m(1) + 2*m(2) + 2*m(3) + m(4))/6.0
 10   continue
      write(13,*) t,x(1),cos(t),cos(t)-x(1)

      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