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


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


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


stop
end

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