Getting NaN in Rung-Kutta 4, Simulation of 10 planets

I am trying to make a model of planets' movement plot it in 3d using Matlab. I used Newton's law with the gravitational force between two objects and I got the differential equation which can be found in the images below.

then I used the numerical method Runge-Kutta-4 to solve it and here's the code:

function [y,t]=RK4(pos,a,b,N)
h=(b-a)/N;
y = zeros(N*10,6);
t = zeros(N,1);
y(1,:)=pos(1,:);
y(2,:)=pos(2,:);
y(3,:)=pos(3,:);
y(4,:)=pos(4,:);
y(5,:)=pos(5,:);
y(6,:)=pos(6,:);
y(7,:)=pos(7,:);
y(8,:)=pos(8,:);
y(9,:)=pos(9,:);
y(10,:)=pos(10,:);
t(1)=a;
for i=1:N
    
    t(i+1)=a+i*h;
    CurrentPos=y((10*i)-9:i*10,:);

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

Where F is the (dY/dt):

function dy=F(t,y,CurrentPos,j)
m=[1.98854E-11 3.302E-18 4.8685E-17 5.97219E-17 6.4185E-18 1.89813E-14 5.68319E-15 8.68103E-16 1.0241E-15 1.307E-19];
G=6.67;
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(i,1)-CurrentPos(j,1));
        deltaY=(CurrentPos(i,2)-CurrentPos(j,2));
        deltaZ=(CurrentPos(i,3)-CurrentPos(j,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

Finally applied the function for the Initial States from JPL HORIZONS System:

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

[yy,t]=RK4(pos,0,100,5);

[James Tursa]

Could you please use comments to indicate the units used for all of your numbers? E.g., why is G listed as 6.67 instead of 6.67e-11? What are the units of the masses and initial positions and velocities?

[Bashar Al-Mohammad]

Thanks for being here,

I use A=10^10m, for the

When using A for distance G became in ordre E-41, Since G and the mass of planet is multiplied by each other I simplified the numbers. Ex: Mass of the Earth: 5.97219E+24 after using A for distance and I multiplied by E-41 then the mass of the Earth: 5.97219E-17.

Thanks again

Answer 1

[James Tursa]

OK, I will take a look at things. But your attempt at "simplifying" your numbers by using a custom distance unit and apparently also custom mass unit only adds confusion IMO. I would advise abandoning this right away and go back to using standard units. In fact, that is the first thing I would do to check this is to convert everything to standard units and run your code as is to see if it is a units problem or a code logic problem.

[Bashar Al-Mohammad]

The new and improve RK4:

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

The new and improved (dY/dt):

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

The new and a bit satisfying results:

>> PlanetaryTest
>> yy

yy =

   1.0e+12 *

    0.0002    0.0010   -0.0000   -0.0000    0.0000    0.0000
   -0.0568   -0.0274    0.0029    0.0000   -0.0000   -0.0000
    0.0428    0.1001   -0.0011   -0.0000    0.0000    0.0000
   -0.1438   -0.0400   -0.0000    0.0000   -0.0000    0.0000
   -0.1147   -0.1963   -0.0013    0.0000   -0.0000   -0.0000
   -0.5669   -0.5775    0.0151    0.0000   -0.0000   -0.0000
    0.0821   -1.5024    0.0229    0.0000    0.0000   -0.0000
    2.6251    1.4027   -0.0288   -0.0000    0.0000    0.0000
    4.3030   -1.2422   -0.0736    0.0000    0.0000   -0.0000
    1.6555   -4.7350    0.0278    0.0000    0.0000   -0.0000
    0.0002    0.0010   -0.0000   -0.0000    0.0000    0.0000
   -0.0568   -0.0274    0.0029    0.0000   -0.0000   -0.0000
    0.0428    0.1001   -0.0011   -0.0000    0.0000    0.0000
   -0.1438   -0.0400   -0.0000    0.0000   -0.0000    0.0000
   -0.1147   -0.1963   -0.0013    0.0000   -0.0000   -0.0000
   -0.5669   -0.5775    0.0151    0.0000   -0.0000   -0.0000
    0.0821   -1.5024    0.0229    0.0000    0.0000   -0.0000
    2.6251    1.4027   -0.0288   -0.0000    0.0000    0.0000
    4.3030   -1.2422   -0.0736    0.0000    0.0000   -0.0000
    1.6555   -4.7350    0.0278    0.0000    0.0000   -0.0000
    0.0002    0.0010   -0.0000   -0.0000    0.0000    0.0000
   -0.0568   -0.0274    0.0029    0.0000   -0.0000   -0.0000
    0.0428    0.1001   -0.0011   -0.0000    0.0000    0.0000
   -0.1438   -0.0400   -0.0000    0.0000   -0.0000    0.0000
   -0.1147   -0.1963   -0.0013    0.0000   -0.0000   -0.0000
   -0.5669   -0.5775    0.0151    0.0000   -0.0000   -0.0000
    0.0821   -1.5024    0.0229    0.0000    0.0000   -0.0000
    2.6251    1.4027   -0.0288   -0.0000    0.0000    0.0000
    4.3030   -1.2422   -0.0736    0.0000    0.0000   -0.0000
    1.6555   -4.7350    0.0278    0.0000    0.0000   -0.0000
    0.0002    0.0010   -0.0000   -0.0000    0.0000    0.0000
   -0.0568   -0.0274    0.0029    0.0000   -0.0000   -0.0000
    0.0428    0.1001   -0.0011   -0.0000    0.0000    0.0000
   -0.1438   -0.0400   -0.0000    0.0000   -0.0000    0.0000
   -0.1147   -0.1963   -0.0013    0.0000   -0.0000   -0.0000
   -0.5669   -0.5775    0.0151    0.0000   -0.0000   -0.0000
    0.0821   -1.5024    0.0229    0.0000    0.0000   -0.0000
    2.6251    1.4027   -0.0288   -0.0000    0.0000    0.0000
    4.3030   -1.2422   -0.0736    0.0000    0.0000   -0.0000
    1.6555   -4.7350    0.0278    0.0000    0.0000   -0.0000

In this time I used (m),(s),(kg). and I got some numbers. I am too tired to plot these but as soon as I do, I will get backk to here and tell. I apperciate your time Mr. Tursa. I aslo appreciate any improvement suggestions on the RK4 function. Thanks again.

[Bashar Al-Mohammad]

Update plotting went wrong. I extracted the x y z coordinates to plot the earth and I got a straight line. 'h' step size used is 1 day variation.