function  dy=kplf(t,r)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
global gamk;
dy=zeros(4,1);
dy(1)=r(3);
dy(2)=r(4);
dy(3)=-r(1)*gamk/sqrt((r(1)^2+r(2)^2))^3;
dy(4)=-r(2)*gamk/sqrt((r(1)^2+r(2)^2))^3;
end

