我试图使行星的运动模型绘制它在3D使用Matlab.我用牛顿定律与两个物体之间的引力，我得到的微分方程如下：
x1c 0d1x的数据
matlab代码：

function dy=F(t,y,CurrentPos,j)
m=[1.98854E+30 3.302E+23 4.8685E+24 5.97219E+24 6.4185E+23 1.89813E+27 5.68319E+26 8.68103E+25 1.0241E+26 1.307E+22];
G=6.67E-11;
dy = zeros(6,1);
dy(1) = y(4);
dy(2) = y(5);
dy(3) = y(6);
for i=1:10
    if i~=j
        deltaX=(CurrentPos(j,1)-CurrentPos(i,1));
        deltaY=(CurrentPos(j,2)-CurrentPos(i,2));
        deltaZ=(CurrentPos(j,3)-CurrentPos(i,3));
        ray=sqrt((deltaX^2)+(deltaY^2)+(deltaZ^2));
        dy(4) = dy(4) + G*m(i)*(deltaX/(ray^3));
        dy(5) = dy(5) + G*m(i)*(deltaY/(ray^3));
        dy(6) = dy(6) + G*m(i)*(deltaZ/(ray^3));
    end
end

字符串
其中m数组是行星的质量。
然后我用数值方法Runge-Kutta-4来解决它，这是代码：

function [y,t]=RK4(F,intPos,a,b,N)

h=(b-a)/N;
t=zeros(N,1);
y = zeros(10*N,6);
y(1,:)=intPos(1,:);
y(2,:)=intPos(2,:);
y(3,:)=intPos(3,:);
y(4,:)=intPos(4,:);
y(5,:)=intPos(5,:);
y(6,:)=intPos(6,:);
y(7,:)=intPos(7,:);
y(8,:)=intPos(8,:);
y(9,:)=intPos(9,:);
y(10,:)=intPos(10,:);

t(1)=a;

for i=1:N
    
    t(i+1)=a+i*h;
    CurrentPos=y((i*10)-9:i*10,:);
%     CurrentPos(1,:)=intPos(1,:);
    y((i*10)+1,:)=intPos(1,:);
    for j=2:10
        k1=F(t(i),y(((i-1)*10)+j,:),CurrentPos,j);
        k2=F(t(i)+h/2,y(((i-1)*10)+j,:)+(h/2).*k1',CurrentPos,j);
        k3=F(t(i)+h/2,y(((i-1)*10)+j,:)+(h/2).*k2',CurrentPos,j);
        k4=F(t(i)+h,y(((i-1)*10)+j,:)+h.*k3',CurrentPos,j);
        y((i*10)+j,:)=y(((i-1)*10)+j,:)+(h/6)*(k1+2*k2+2*k3+k4)';
    end
end

型
最后应用JPL HORIZONS系统的初始状态函数：

format short

intPos=zeros(10,6);
intPos(1,:)=[1.81899E+08 9.83630E+08 -1.58778E+07 -1.12474E+01 7.54876E+00 2.68723E-01];
intPos(2,:)=[-5.67576E+10 -2.73592E+10 2.89173E+09 1.16497E+04 -4.14793E+04 -4.45952E+03];
intPos(3,:)=[4.28480E+10 1.00073E+11 -1.11872E+09 -3.22930E+04 1.36960E+04 2.05091E+03];
intPos(4,:)=[-1.43778E+11 -4.00067E+10 -1.38875E+07 7.65151E+03 -2.87514E+04 2.08354E+00];
intPos(5,:)=[-1.14746E+11 -1.96294E+11 -1.32908E+09 2.18369E+04 -1.01132E+04 -7.47957E+02];
intPos(6,:)=[-5.66899E+11 -5.77495E+11 1.50755E+10 9.16793E+03 -8.53244E+03 -1.69767E+02];
intPos(7,:)=[8.20513E+10 -1.50241E+12 2.28565E+10 9.11312E+03 4.96372E+02 -3.71643E+02];
intPos(8,:)=[2.62506E+12 1.40273E+12 -2.87982E+10 -3.25937E+03 5.68878E+03 6.32569E+01];
intPos(9,:)=[4.30300E+12 -1.24223E+12 -7.35857E+10 1.47132E+03 5.25363E+03 -1.42701E+02];
intPos(10,:)=[1.65554E+12 -4.73503E+12 2.77962E+10 5.24541E+03 6.38510E+02 -1.60709E+03];

[yy,t]=RK4(@F,intPos,0,1e8,1e3);
x=zeros(101,1);
y=zeros(101,1);
z=zeros(101,1);
for i=1:1e3
    x(i,:)=yy((i-1)*10+4,1);
    y(i,:)=yy((i-1)*10+4,2);
    z(i,:)=yy((i-1)*10+4,3);
end

plot3(x,y,z)

型
最后，结果一点也不令人满意，我得到了很多“NAN”，然后我对RK 4方法做了一些调整，开始得到数字，但当我绘制它们时，我绘制的是一条线而不是轨道。
任何帮助都将不胜感激。先谢了。