real y(2),work(15),m,k
integer iwork(5)
common omega2,beta
external func
open(10,file="harm.dat")
k=5.
m=0.1
omega2=k/m
alpha=0.1
beta=alpha/(2.*m)
relerr=1.e-8
abserr=1.e-8
neqn=2
t=0.
iflag=1
y(1)=0.1
y(2)=0.
nt=10000
dt=0.001
too=0.
neqn=2
x=y(1)

do i=1,nt
tout=t+dt
call rkf45(func,neqn,y,t,tout,relerr,abserr,iflag,work,iwork)
if (iflag.ne.2) then
iflag=2
endif

dx=y(1)-x
too=too-alpha*y(2)*dx
x=y(1)
ekin=0.5*m*y(2)**2
epot=0.5*k*y(1)**2
etot=ekin+epot
ecorr=etot-too

write(10,'(100f15.5)') tout,y,ekin,epot,etot,ecorr
enddo


stop
end

subroutine func(t,y,dy)
common omega2,beta
real y(2),dy(2)
dy(1)=y(2)
dy(2)=-omega2*y(1)-2.*beta*y(2)
return
end