function  dy=kplf2(t,y)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
global gm gk omegak kor fi0;
dy=zeros(4,1);
xk=kor*cos(omegak*t+fi0);
yk=kor*sin(omegak*t+fi0);
xx=(y(1)-xk);
yy=(y(2)-yk);
rr=sqrt(xx^2+yy^2);
dy(1)=y(3);
dy(2)=y(4);
dy(3)=-y(1)*gm/sqrt((y(1)^2+y(2)^2))^3-xx*gk/rr^3;
dy(4)=-y(2)*gm/sqrt((y(1)^2+y(2)^2))^3-yy*gk/rr^3;
end

