program enshu14rkg2f                                 
implicitnone
integer,parameter::kosuu=16
real,dimension(0:kosuu)::t,vrkg,xrkg
real::dt,kv,kx,qv,qx
real,parameter::t_kaishi=0.0,t_shuryo=2.0
integer::i
real,parameter::pi=4.0*atan(1.0),pi2jo=pi**2
dt=(t_shuryo-t_kaishi)/real(kosuu)
doi=0,kosuu
t(i)=t_kaishi+dt*real(i)
enddo
vrkg(0)=pi
xrkg(0)=0.0
doi=0,kosuu-1
kv=dt*dvdt(xrkg(i))
kx=dt*dxdt(vrkg(i))
vrkg(i+1)=vrkg(i)+0.5*kv
xrkg(i+1)=xrkg(i)+0.5*kx
qv=kv
qx=kx
kv=dt*dvdt(xrkg(i+1))
kx=dt*dxdt(vrkg(i+1))
vrkg(i+1)=vrkg(i+1)+(1.0-sqrt(0.5))*(kv-qv)
xrkg(i+1)=xrkg(i+1)+(1.0-sqrt(0.5))*(kx-qx)
qv=(3.0*sqrt(0.5)-2.0)*qv+(2.0-sqrt(2.0))*kv
qx=(3.0*sqrt(0.5)-2.0)*qx+(2.0-sqrt(2.0))*kx
kv=dt*dvdt(xrkg(i+1))
kx=dt*dxdt(vrkg(i+1))
vrkg(i+1)=vrkg(i+1)+(1.0+sqrt(0.5))*(kv-qv)
xrkg(i+1)=xrkg(i+1)+(1.0+sqrt(0.5))*(kx-qx)
qv=-(2.0+3.0*sqrt(0.5))*qv+(2.0+sqrt(2.0))*kv
qx=-(2.0+3.0*sqrt(0.5))*qx+(2.0+sqrt(2.0))*kx
kv=dt*dvdt(xrkg(i+1))
kx=dt*dxdt(vrkg(i+1))
vrkg(i+1)=vrkg(i+1)+(0.5*kv-qv)/3.0
xrkg(i+1)=xrkg(i+1)+(0.5*kx-qx)/3.0
enddo
doi=0,kosuu
write(*,*)t(i),vrkg(i),xrkg(i)
enddo
contains
functiondvdt(x)
real,intent(in)::x
real::dvdt
dvdt=-pi2jo*x
endfunctiondvdt
functiondxdt(v)
real,intent(in)::v
real::dxdt
dxdt=v
endfunctiondxdt
endprogramenshu14rkg2f