gravity simulation in matlab

main.m
% Parameters
num_bodies = 3;
G = 6.6743e-11; % gravitational constant
dt = 0.01; % time step

% Initial conditions
r0 = [1, 0, 0; -1/2, sqrt(3)/2, 0; -1/2, -sqrt(3)/2, 0]; % initial positions
v0 = [0, 0, 0; 0, 0, 0; 0, 0, 0]; % initial velocities
m = [5.972e24; 7.342e22; 1]; % masses

% Simulation loop
r = r0; v = v0; % set initial conditions
for i = 1:1000 % run for 1000 time steps
    % Compute acceleration using current positions
    a = zeros(num_bodies, 3);
    for j = 1:num_bodies
        for k = 1:num_bodies
            if j ~= k
                rjk = r(j,:) - r(k,:);
                a(j,:) = a(j,:) - G*m(k)*rjk/norm(rjk)^3;
            end
        end
    end
    
    % Update positions using Verlet algorithm
    r_new = 2*r - r0 + a*dt^2;
    
    % Update velocities using central difference method
    v_new = (r_new - r0)/(2*dt);
    
    % Update variables for next iteration
    r0 = r; r = r_new;
    v0 = v; v = v_new;
    
    % Plot current positions
    plot3(r(:,1), r(:,2), r(:,3), 'o');
    axis([-2, 2, -2, 2, -2, 2]*10^8);
    drawnow;
end
1074 chars
40 lines