Program harmony
c-----------------------------------------------------------------
c October 2006
c Harmonic oscillator (can be modified for any 1D second order ODE
c Program may call various ODU solvers
c-----------------------------------------------------------------
	Real*8 t0, x0, v0, ei
	Real*8 d1x,d2x,ti,xi,vi,tf,xf,vf
	Real*8 tmax, k, m
	Integer*4 i
	Common/const/k,m
	External d1x, d2x
      open (unit=7,file="harmony0.dat")
c* parameters for oscillator
      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.02
	tmax = 20.00

	ti = t0
      xi = x0
	vi = v0
	ei = vi**2/2 + xi**2/2
	do while(ti.le.tmax)
	   tf = ti + dt

c*** select ONE of the methods below

c*** Runge-Kutta
c	   call rk4_2d(d1x,d2x,ti,xi,vi,tf,xf,vf)
c*** Euler: predictor-corrector
c	   call euler2m(d1x,d2x,ti,xi,vi,tf,xf,vf)
c*** simple Euler
           call euler2d(d1x,d2x,ti,xi,vi,tf,xf,vf)
	   energy = vf**2/2 + xf**2/2
         write (7,100) tf,xf,vf, energy
	   ti = tf
	   xi = xf
	   vi = vf
	end do
100	format(4(1pe12.4))
	stop
	End

	Function d1x(t,x,v)
c--------------------------------------
c function dx/dt
c--------------------------------------
	Real*8 d1x, t, x, v
	  d1x = v
	return
	end

	Function d2x(t,x,v)
c--------------------------------------
c function d2x/dt2
c--------------------------------------
	Real*8 d2x, t, x, v, k, m
	Common/const/k,m
        d2x = (-1.0*k/m)*x
	return
        end