一道matlab题目，求解答。

Question

Answer 1

你好，可以参考为你写的代码（如有帮助还望题主给个宝贵的采纳支持一下答主哟）：
微分方程：

Matlab代码（有动图）

% x(1)对应r，x(2)对应dr/dt，x(3)对应theta，x(4)对应dtheta/dt
M=5.965e24; %地球质量
G=6.67e-11; %万有引力常数
odefun = @(t,x)[
    x(2);
    x(1)*x(4)^2-G*M/x(1)^2;
    x(4);
    -2*x(2)*x(4)/x(1);
    ];
r0 = 384403.9e3; %地月距离
dr0 = 0;%初始径速度
theta0 = 0;%初始角度
dtheta0 = 1020/r0;%初始角速度
x0 = [r0,dr0,theta0,dtheta0];
[t,y] = ode45(odefun,[0,3600*24*30*3],x0);%模拟三个月
r = y(:,1);
theta = y(:,3);
xp = r.*cos(theta);
yp = r.*sin(theta);
% 画图
for i = 1:1:numel(xp)
figure(1);clf
plot(xp(1:i)/1e7,yp(1:i)/1e7,'b-')
hold on;
plot(0,0,'bo','markersize',15,'markerfacecolor','b')
plot(xp(i)/1e7,yp(i)/1e7,'mo','markersize',5,'markerfacecolor','m')
axis([-1,1,-1,1]*1.2*r0/1e7)
xlabel('X(万公里)');ylabel('Y(万公里)')
text(xp(i)/1e7+3,yp(i)/1e7+3,'月球')
text(3,3,'地球')
title(['第',num2str(t(i)/3600/24),'天'])
axis equal
if(mod(i,3600)==1)
    pause(0.00001)
end
end

图片效果：