real y(2) , work(15),mass,k integer iwork(5) common beeta, omega0 external f data pi/3.1415926536/ open (10,file='sumb.in') open (20, file='sumb.out') neqn=2 read(10,*)relerr,abserr,beeta,mass,k,dt,nt,y omega0=k/mass talg=0 tlopp=0 iflag=1 call rkf45(f,neqn,y,talg,tlopp,relerr,abserr,iflag,work,iwork) do i =1,nt tlopp=talg+dt iflag=2 call rkf45(f,neqn,y,talg,tlopp,relerr,abserr,iflag,work,iwork) talg=tlopp write(20,100) tlopp,y,mass*y(2)**2/2.,1.-k*cos(2.*pi*y(1)) * ,mass*y(2)**2/2.+1.-k*cos(2.*pi*y(1)) 100 format(6f15.5) enddo stop end subroutine f(t,y,dy) real y(2),dy(2) common beeta, omega0 data pi/3.1415926536/ dy(1)=y(2) dy(2)=-beeta*y(2)-2.*pi*omega0*sin(2.*pi*y(1)) return end