implicit none
real y(2),work(15),l,g,alpha,relerr,abserr,m
real omega2,beta,t,tout,dt,too,ekin,epot,etot
integer iwork(5),iflag,neqn,i,nt
common omega2,alpha
external func
open(10,file="harm_math.dat")
l=5.
g=9.814
omega2=g/l
alpha=0.1
relerr=1.e-8
abserr=1.e-8
neqn=2
t=0.
iflag=1
y(1)=0.
y(2)=4.
nt=10000
dt=0.01
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
enddo


stop
end

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