Program harmony2
c-----------------------------------------------------------------
c Driver program to rkf45.f
c Solution of second-order ODE (initial value problem)
c Application to a simple harmonic oscillator
c (can be easily modified for any 1D second order ODE
c-----------------------------------------------------------------
      Real*8 t0, x0, v0, ei
      Real*8 tmax, k, m
c* for Rkf45
      Real*8 t, x(2), tout, relerr, abserr
      Real*8 work(15)
      integer*4 iwork(5), iflag, neqn
      parameter (relerr=1.0e-9, abserr=0.0)
c* end Rkf45
      Integer*4 i
      Common/const/k,m
      External hoscc
      open (unit=7,file="harmony4.dat")
c* parameters for oscillator
      neqn = 2
      k=1.0
      m=1.0
c* initial conditions
      t0 = 0.0
      x0 = 0.0
      v0 = 1.0
c* time step and max time
      dt   =  0.2
      tmax = 20.0

      t = t0
      x(1) = x0
      x(2) = v0
      iflag = 1
      tout = t
      ei = v0**2/2 + x0**2/2
      do while(tout.le.tmax+dt)
c*** Rkf45
       call rkf45(hoscc,neqn,x,t,tout,relerr,abserr,iflag,work,iwork)
c*********
       ei = x(1)**2/2 + x(2)**2/2
       write (7,100) tout, x(1),x(2), ei
       tout = t+dt
      end do
100   format(5(1pe12.4))
      stop
      end

      Subroutine hoscc(t,x,dx)
c-------------------------------------------------
c definition of the differential equations
c x(1) = x    dx(1)=v (or x(2))
c x(2) = v    dx(2)=dv/dt = -(k/m)*x = -(k/m)*x(1)
c-------------------------------------------------
      Real*8 t, x(2), dx(2), k, m
      Common/const/k,m
      dx(1) = x(2)
      dx(2) = (-1.0*k/m)*x(1)
      return
      end