real y(2),work(15),l,g,m
integer iwork(5)
common omega2,beta
external func
open(10,file="harm.dat")
l=0.1
g=9.814
omega2=l/g
m=0.1
alpha=0.001
beta=alpha/m

relerr=1.e-8
abserr=1.e-8
neqn=2
t=0.
iflag=1
y(1)=0
y(2)=0.3
nt=100000
dt=0.01
tout=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
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))-beta*y(2)
return
end