function du = satellite (t, u, GM)
r =  sqrt(u(1)^2+u(3)^2);
du(1,1) = u(2); du(2,1)=-GM/(r^3)*u(1);
du(3,1) = u(4); du(4,1)=-GM/(r^3)*u(3);
endfunction
