clear all;close all;G=6.67*10-11;% The universal gravitational constantm = 1.989e30,3.5844e23,4.89868e24,5.974e24,6.5714e23,1.89854e27,5.68725e26,8.72204e25,1.02753e26;% An array of massesn = length(m); Number_of_planets = n% The number of massesx_p = 0,58340100000,1.07705e11,1.4959e11,2.27377e11,7.77868e11,1.42709e12,2.87512e12,4.49668e12;% An array of x positionsy_p= 0,0,0,0,0,0,0,0,0;% An array of y positionsx_v = 0,0,0,0,0,0,0,0,0;% An array of x velocitiesy_v = 0,47856.46,34961.72,29780,23883.56,12924.52,9618.94,6789.84,5419.96;% An array of y velocitiesT= 365*24*60*60;%(2*pi*G*m(1)/(x_v(1)3);dt=360*60*60*10;colordef black;ph = plot(x_p,y_p,.,MarkerSize,30);xlabel(distance(km);ylabel(distance(km);title(Planetary motion);for a=0:dt:1000*T%loop for all masses with respect to timefor k=1:n%loop for individual masses calculating neccessary quantitiesdx = x_p - x_p(k);% difference in x positionsdy = y_p - y_p(k);% difference in y positionsmag = (dx.2 + dy.2).0.5;% magnitude of the distance between the 2 massesf = (G*m(k)*m)./(mag.2);% total forcefx_old =(f.*dx./mag);fx_old(k)= 0;fx = sum (fx_old);%Summing the total force in x direction for each planetfy_old =(f.*dy./mag);fy_old(k)= 0;fy = sum (fy_old);%Summing the total force in y direction for each planeta_x = fx/m(k);% calculating acceleration change in x directiona_y = fy/m(k);% calculating acceleration change in y directionx_v(k) = x_v(k) + dt*a_x;% calculating velocity change in x directiony_v(k) = y_v(k) + dt*a_y;% calculating velocity change in y directionx_p_new(k)= x_p(k)+dt*x_v(k);% calculating new x positionsy_p_new(k)= y_p(k)+dt*y_v(k);% calculating new y positionsendx_p=x_p_new;y_p=y_p_new;% overwriting old positions with new positionspause(0.1)plot(x_p,y_p,.,MarkerSize,30)axis(-(1015) 1015 -(1015) 1015)xlabel(distance(km);ylabel(distance(km);title(Planetary motion);drawnow% Plotting graphend