caculate the path for a rocket between earth and moon in matlab

main.m
% Define constants
G = 6.67430e-11; % gravitational constant
Me = 5.972e24; % mass of Earth
Mm = 7.342e22; % mass of Moon
Re = 6.3781e6; % radius of Earth
Rm = 1.7371e6; % radius of Moon
De = 3.844e8; % distance from Earth to Moon

% Define initial conditions
m = 10000; % mass of rocket (kg)
ve = 2000; % initial velocity of rocket (m/s)
xe = Re; % initial position of rocket (m)

% Define thrust properties
thrust = 5000; % thrust of rocket (N)
isp = 300; % specific impulse of rocket (s)

% Define equations of motion
f = @(t, y) [
    y(4); % x velocity
    y(5); % y velocity
    y(6); % z velocity
    -G*Me*y(1)/(norm(y(1:3))^3) - G*Mm*(y(1:3)-[De;0;0])/(norm(y(1:3)-[De;0;0])^3) + thrust*y(7)*y(4)/m; % x acceleration
    -G*Me*y(2)/(norm(y(1:3))^3) - G*Mm*y(2:4)/(norm(y(2:4))^3) + thrust*y(7)*y(5)/m; % y acceleration
    -G*Me*y(3)/(norm(y(1:3))^3) - G*Mm*y(3:5)/(norm(y(3:5))^3) + thrust*y(7)*y(6)/m; % z acceleration
    -thrust/(m*isp)]; % fuel consumption rate

% Integrate equations of motion
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-6);
[t, y] = ode45(f, [0, 3500000], [xe; 0; 0; 0; ve; 0; 1], options);

% Plot path of rocket
plot3(y(:,1), y(:,2), y(:,3), 'linewidth', 2)
hold on
tmoon = linspace(0, 2*pi, 100);
xmoon = Rm*cos(tmoon) + De;
ymoon = Rm*sin(tmoon);
plot3(xmoon, ymoon, zeros(size(xmoon)), 'k', 'linewidth', 2)
hold off
xlabel('x (m)')
ylabel('y (m)')
zlabel('z (m)')
axis equal
1415 chars
44 lines