real y(4),mass,work(27),k
integer iwork(5)
common param,r0
external func
open(10,file="ruther.dat")
b=0.05
r0=1.
q1=3.e-6
q2=0.1e-6
k=9.e9
mass=1.e-6

param=q1*q2*k/mass
relerr=1.e-5
abserr=1.e-5
dt=0.0001
nt=100
iflag=1

y(1)=0.
y(2)=b
y(3)=50.
y(4)=0.
neqn=4
t=0.


do it=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,*) tout,y
enddo

stop
end

subroutine func(t,y,yp)
real y(4),yp(4)
common param,r0

yp(1)=y(3)
yp(2)=y(4)
yp(3)=param*(y(1)-r0)/(sqrt((y(1)-r0)**2+y(2)**2))**3
yp(4)=param*y(2)/(sqrt((y(1)-r0)**2+y(2)**2))**3
return
end