implicit none
real y(2),work(15),relerr,abserr,m,alpha
real A,beta,t,tout,dt,ekin,epot,etot,gamma
real u0,pi,L
integer iwork(5),iflag,neqn,i,nt
common A,beta,gamma
external func
open(10,file="harm_math.dat")
pi=3.1415926536
alpha=0.1
m=1.
L=1.
u0=1.

A=u0*pi/(m*L)
beta=alpha/m
gamma=2*pi/L

relerr=1.e-9
abserr=1.e-9
neqn=2
t=0.
iflag=1
y(1)=0.
y(2)=2
nt=10000
dt=0.01
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 A,beta,gamma
real y(2),dy(2)
dy(1)=y(2)
dy(2)=-A*sin(gamma*y(1))-beta*y(2)
return
end