模拟并显示太阳系的形成

模拟并显示太阳系的形成。
① 采用“粒子”作为基本计算单元，基本粒子在距离足够近 时结合成大的物体。
当质量增加到不同级别后发不同颜色 的光（可自行设定颜色变化图）；
粒子数量多时，采用并 行计算加速
② 物理模型：仅利用万有引力公式𝑭 = 𝑮 𝑴𝒎 𝒓 𝟐
③ 提示：；仿真中时间间隔自定。

你好，可以尝试我为你写的代码：

clc;clear
rng(1)
n = 400;
m = ones(n,1)*10e24; % 初始质量kg
G = 6.754e-11; % 引力常数N·m²/kg²
x = rand(n,1)*1e14; % 初始坐标m
y = rand(n,1)*1e14; % 初始坐标m
vx = zeros(n,1);
vy = zeros(n,1);
dt = 5e8; % 时间步
for tt = 1:1000000
    n = numel(x);
    fx = zeros(n,1);
    fy = zeros(n,1);
    del = [];
    for i = 1:n
        for j = i+1:n
            dx = x(j) - x(i);
            dy = y(j) - y(i);
            m2 = m(i)*m(j);
            r2 = dx*dx + dy*dy;
            r = sqrt(r2);
            if(r2<(1e12)^2 && m(i)~=0 && m(j)~=0)
                vx(i) = (m(i)*vx(i)+m(j)*vx(j))/(m(i)+m(j));
                vy(i) = (m(i)*vy(i)+m(j)*vy(j))/(m(i)+m(j));
                vx(j) = 0;
                vy(j) = 0;
                m(i) = m(i)+m(j);
                m(j) = 0;
                del = [del, j];
            else
                fx(i) = fx(i) + G*m2/r2*dx/r;
                fy(i) = fy(i) + G*m2/r2*dy/r;
                fx(j) = fx(j) - G*m2/r2*dx/r;
                fy(j) = fy(j) - G*m2/r2*dy/r;
            end
        end
        
    end
    % 删除被吞并的星体
    m(del) = [];
    fx(del) = [];
    fy(del) = [];
    vx(del) = [];
    vy(del) = [];
    x(del) = [];
    y(del) = [];
    ax = fx./m;
    ay = fy./m;
    ax(isnan(ax)) = 0;
    ay(isnan(ay)) = 0;
    vx = vx + ax*dt;
    vy = vy + ay*dt;
    x = x + vx*dt;
    y = y + vy*dt;
    if(mod(tt,10)==0)
        figure(1);clf
        scatter(x,y,10*(m/5e24).^(1/2),m, 'filled')
        colormap(jet(8))
        h = colorbar;
        set(get(h,'Title'),'string','kg');
        axis([0,1,0,1]*1e14)
        axis equal
        grid on
        title([num2str(dt*tt/(365*3600*24)),'年'])
        pause(0.00000001)
    end
end

融合效果

可以设计成三维立体的吗？并且增大粒子的数目？其他部分都可以说是完美，尤其是可以改变大小和颜色，更加的优美与直观。如果可以解决前两个问题就更好了，真的是很感谢。

好的，谢谢您抽空帮我完成这些，万分感谢

很抱歉，以我的能力只能将物理关系改成三维，但是无法运用MATLAB将粒子运动过程通过三维的方式，三维的图像表达出来，很抱歉，能帮助修改一下算法吗？
我也会努力学习的。