function f=rk(t,w)
x=w(1);
xdot=w(2);
y=w(3);
ydot=w(4);
xp=t/5;
r=sqrt((x-xp)^2+y^2);
force=1/r^2;
ct=(x-xp)/r;
st=y/r;
f=[xdot; -force*ct; ydot; -force*st];

