dotu = odefunc(t,u)
    x = u(1:N)
    y = u(N+1,2*N)
    ....

    r2 = x.^2+y.^2
    r = r2.^0.5
    ....
    doty = gamma*(mu-r2).*y+omega.*x
    Ft=F(t)
    dotomega = -epsilon*Ft*y./r
    ...