program enshu14rkg2                                  
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*pi2jo*xrkg(i)
kx=dt*vrkg(i)
vrkg(i+1)=vrkg(i)+0.5*kv
xrkg(i+1)=xrkg(i)+0.5*kx
qv=kv
qx=kx
kv=-dt*pi2jo*xrkg(i+1)
kx=dt*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*pi2jo*xrkg(i+1)
kx=dt*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*pi2jo*xrkg(i+1)
kx=dt*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
endprogramenshu14rkg2