% Моделирование вращения Земли вокруг Солнца
function earth
    rz= 149600000*1000; % Радиус (начальное расстояние)
    vz= 29783;  % Начальная скорость
    options= odeset('AbsTol',1e+1,...
        'RelTol',1e-4); % Здесь для некоторых версий матлаба
                        % нужна настройка. Иначе виснет
    [t,xy]= ode45(@erp,...
        [0 365.25*24*60*60],... % интервал инт-я - 1 год
        [rz 0 0 vz],options); % В начальный момент Земля на оси x
    plot(xy(:,1),xy(:,3));
    figure
    plot(t,xy(:,1),t,xy(:,3));
    figure
    plot(t,xy(:,2),t,xy(:,4)...
        ,t,sqrt(xy(:,2).^2+xy(:,4).^2));

end
function dxy=erp(~,z)
    x= z(1); vx= z(2); 
    y= z(3); vy= z(4);
    e= [-x -y]./sqrt(x^2+y^2); % орт вектора, напр. из Земли на Солнце
    G= 6.67e-11;
    ms= 1.989e+30;  % масса Солнца
    a= e*G*ms./(x^2+y^2); % вектор ускорения
    dxy= [vx;a(1);vy;a(2)];
end