simulate the solar system orbits in matlab

main.m
% define constants
G = 6.67430e-11;    % gravitational constant
M_sun = 1.989e30;   % mass of the sun

% define initial conditions (position and velocity) for each planet
% (in meters and meters/second, respectively)
S = [0; 0; 0];
V = [0; 0; 0];
M = 5.972e24;
E = [1.49597871e11; 0; 0];
V_E = [0; 29.29e3; 0];
M_sun = 1.9891e30;
M_mars = 6.41693e23;
R_mars = [227.936e9; 0; 0]; 
V_mars = [0; 24.077e3; 0];

% define time step and duration of simulation
dt = 3600*24; % one day
t_end = 3600*24*365.25*10; % ten years

% initialize arrays to store position and velocity data for each planet
num_steps = t_end/dt;
pos_S = zeros(3,num_steps+1); vel_S = zeros(3,num_steps+1);
pos_E = zeros(3,num_steps+1); vel_E = zeros(3,num_steps+1);
pos_Mars = zeros(3,num_steps+1); vel_Mars = zeros(3,num_steps+1);

% loop over time steps and update position and velocity of each planet
for i = 1:num_steps+1
    % update positions of planets
    pos_S(:,i) = S; pos_E(:,i) = E; pos_Mars(:,i) = R_mars;
    
    % calculate acceleration of each planet due to gravitational forces
    r_SE = E - S;
    r_SMars = R_mars - S;
    r_EMars = R_mars - E;
    a_S = G*M_sun/norm(S)^3 * S + G*M_mars/norm(r_SMars)^3 * r_SMars;
    a_E = G*M_sun/norm(r_SE)^3 * r_SE + G*M_mars/norm(r_EMars)^3 * r_EMars;
    a_Mars = G*M_sun/norm(R_mars)^3 * R_mars + G*M/norm(r_SMars)^3 * (-r_SMars) + G*M_e/norm(r_EMars)^3 * (-r_EMars);
    
    % update velocities of planets
    vel_S(:,i) = V + a_S*dt;
    vel_E(:,i) = V_E + a_E*dt;
    vel_Mars(:,i) = V_mars + a_Mars*dt;
    
    % update positions of planets
    S = S + vel_S(:,i)*dt;
    E = E + vel_E(:,i)*dt;
    R_mars = R_mars + vel_Mars(:,i)*dt;
end

% plot trajectories of each planet
figure; hold on;
plot3(pos_S(1,:),pos_S(2,:),pos_S(3,:),'b');
plot3(pos_E(1,:),pos_E(2,:),pos_E(3,:),'g');
plot3(pos_Mars(1,:),pos_Mars(2,:),pos_Mars(3,:),'r');
xlabel('x (m)'); ylabel('y (m)'); zlabel('z (m)');
legend('Sun','Earth','Mars');
1952 chars
58 lines