Planet Orbit- Dormand Prince Algorithm

function PlanetOrbit
[S,M,AU,C1,C2,R]=StateInitialization;

dt=3600;
t=0;
DT=datenum(2019,12,1,0,0,0);
SF=1.157407407407407e-05; %keeps track of date
plotState(S,AU,C1,C2,R,DT);

eps=0.000001; %error allowance in one step calculation.
    t0=0; %initiallization of variables.
    y0=0;
    h0=0.1;
    while(t0<=tf) %tf is the ending time of the calculation.
    
      k1=h*f(t0,y0);
      k2=h*f(t0 + (1/5)*h, y0 + (1/5)*k1);
      k3=h*f(t0 + (3/10)*h, y0+(3/40)*k1+(9/40)*k4);
      k4=h*f(t0+(4/5)*h, y0 + (44/45)*k1-(56/12)*k2+(32/9)*k3);
      k5=h*f(t0+(8/9)*k1, y0+(19372/6561)*k1-(25360/2187)*k2+(64448/6561)*k3-(212/729)*k4);
      k6=h*f(t0+h, y0+(9017/3168)*k1-(355/33)*k2-(46732/5247)*k3+(49/176)*k4-(5103/18656)*k5);
      k7=h*f(t0+h,y0+(35/384)*k1+(500/1113)*k3+(125/192)*k4-(2187/6784)*k5+(11/84)*k6);
      y1 = y0 + 35/384*k1 + 500/113*k3 + 125/192*k4 - 2187/6784*k5 +11/84*k6;
      z1 = y0 + 5179/57600*k1 + (7571/16695)*k3+(393/640)*k4-(92097/339200)*k5+(187/2100)*k6+(1/40)*k7; %to estimate the error.
      err = abs(z1-y1); %error estimation
      s=pow(eps*h0/(2*err),1/5);
      h1=s*h0; %optimal time interval. used in the next step.
      if(h1<hmin)
          h1=hmin; %confine time interval between hmin and hmax.
      else if(h1>hmax)h1=hmax;
      t0=t0+h0; 
      y0=y1;
      h0=h1;
          end
      end
    end
    

for j=1:50000
    S=S+dt*physics(S,M);
    t=t+dt;
    DT=DT+SF*dt;
    if mod(j,10)==0
        plotState(S,AU,C1,C2,R,DT);
    end
end
end

function plotState(S,AU,C1,C2,R,DT)
    NB=round(length(S)/6);
    clf;
    hold on;
    for j=1:NB
        bi=6*(j-1)+1;
        plot3(S(bi)/AU,S(bi+1)/AU,S(bi+2)/AU,C1(j,:),'MarkerSize',R(j),'MarkerFaceColor',C2(j));
    end
    axis([-2,2,-2,2,-2,2])
    title(datestr(round(DT,4)))
    drawnow;
end

I get the error message:

>> PlanetOrbit

Index exceeds matrix dimensions.

Error in PlanetOrbit>plotState (line 59)

plot3(S(bi)/AU,S(bi+1)/AU,S(bi+2)/AU,C1(j,:),'MarkerSize',R(j),'MarkerFaceColor',C2(j));

Error in PlanetOrbit (line 12)

plotState(S,AU,C1,C2,R,DT);

The other functions used are:

%start with this and finish this
function [S,M,AU,C1,C2,R]=StateInitialization

% S state of system (x,y,z,vx,vy,vz) meters and m/s long vector
% M list of masses (sun,earth, etc.)
% AU astronomical unit(meters) earth to sun
% C1 keeps track of outer circle, fill interior blue
% C2 outer color
% R how big dots appear

AU = 149597870700; 

C1(1,:)='yo'; C2(1)='y';
C1(2,:)='bo'; C2(2)='b';
C1(3,:)='mo'; C2(3)='m';
C1(4,:)='go'; C2(4)='g';
C1(5,:)='ro'; C2(5)='r';
C1(6,:)='ko'; C2(6)='k';
C1(7,:)='co'; C2(7)='c';
C1(8,:)='wo'; C2(8)='w';

R=[12;8;8;5];

S=zeros(6,1);
%Sun
X=0; Y=0; Z=0;
VX=0; VY=0; VZ=0;
j=1; bi=6*(j-1)+1;
S(bi)=1000*X; S(bi+1)=1000*Y; S(bi+2)=1000*Z;
S(bi+3)=1000*VX; S(bi+4)=1000*VY; S(bi+5)=1000*VZ;
M(1)=1.989*1e30;

%Earth/Moon
 X = 5.476756046031971E+07; Y = 1.369854390984182E+08; Z =-6.418494023710489E+03;
 VX=-2.814628809603848E+01; VY= 1.094601049527264E+01; VZ=-4.343632416126120E-04;
j=j+1; bi=6*(j-1)+1;
S(bi)=1000*X; S(bi+1)=1000*Y; S(bi+2)=1000*Z;
S(bi+3)=1000*VX; S(bi+4)=1000*VY; S(bi+5)=1000*VZ;
M(2)=5.972*1e24+7.34767309*1e22;

%Venus
X = 7.690918093481298E+07; Y =-7.695356068929358E+07; Z =-5.494117392182905E+06;
VX= 2.454395327263311E+01; VY= 2.462723606903867E+01; VZ=-1.078441904000677E+00;

j=j+2; bi=6*(j-1)+1;
S(bi)=1000*X; S(bi+1)=1000*Y; S(bi+2)=1000*Z;
S(bi+3)=1000*VX; S(bi+4)=1000*VY; S(bi+5)=1000*VZ;
M(3)=4.867*1e24;

%Apophis
X = 1.006188535413713E+08; Y = 7.025826519631593E+07; Z =-1.349533159624398E+06;
VX=-1.511473674904992E+01; VY= 3.112188699288047E+01; VZ=-2.016175343620677E+00;
j=j+3; bi=6*(j-1)+1;
S(bi)=1000*X; S(bi+1)=1000*Y; S(bi+2)=1000*Z;
S(bi+3)=1000*VX; S(bi+4)=1000*VY; S(bi+5)=1000*VZ;
M(4)=26.99*1e9;

%Saturn
X = 5.454260468832656E+08; Y = -1.399004902083351E+09; Z = 2.611933568101764E+06;
VX= 8.479463517851789E+00; VY= 3.483787999150237E+00; VZ=-3.985176203231436E-01;
j=j+4; bi=6*(j-1)+1;
S(bi)=1000*X; S(bi+1)=1000*Y; S(bi+2)=1000*Z;
S(bi+3)=1000*VX; S(bi+4)=1000*VY; S(bi+5)=1000*VZ;
M(5) = 5.6834*1e26;

%Jupiter
X = 4.421300143790728E+07; Y =-7.824664674728345E+08; Z = 2.260798916549563E+06;
VX= 1.290169373776935E+01; VY= 1.355135395534917E+00; VZ=-2.942549827353074E-01;
j=j+5; bi=6*(j-1)+1;
S(bi)=1000*X; S(bi+1)=1000*Y; S(bi+2)=1000*Z;
S(bi+3)=1000*VX; S(bi+4)=1000*VY; S(bi+5)=1000*VZ;
M(6) = 1.89813*1e27;

%Neptune
X = 4.371602845443687E+09; Y =-9.668331843241746E+08; Z =-8.085109646501470E+07;
VX= 1.151818638737934E+00; VY= 5.342596132371595E+00; VZ=-1.370365801570574E-01;
j=j+6; bi=6*(j-1)+1;
S(bi)=1000*X; S(bi+1)=1000*Y; S(bi+2)=1000*Z;
S(bi+3)=1000*VX; S(bi+4)=1000*VY; S(bi+5)=1000*VZ;
M(7) = 1.02413*1e26;

%Uranus
X = 2.437829683225929E+09; Y = 1.688160852693148E+09; Z =-2.530390129291379E+07;
VX=-3.914948553525651E+00; VY= 5.282609057004567E+00; VZ= 7.012996969005503E-02;
j=j+7; bi=6*(j-1)+1;
S(bi)=1000*X; S(bi+1)=1000*Y; S(bi+2)=1000*Z;
S(bi+3)=1000*VX; S(bi+4)=1000*VY; S(bi+5)=1000*VZ;
M(8) = 8.6813*1e25;
 
end

and

function dzdt = physics(z,t)
    dzdt=0*z;
    dzdt(3) = z(1) * cos(z(2));
    dzdt(4) = z(1) * sin(z(2));
    dzdt(2) = -9.81/z(1) * cos(z(2));
    D = ((0.5) * (1.29)* (0.72))/2 * (z(1)^2);
    dzdt(1) = -D/(80) - 9.81*sin(z(2));
end

function[T] = GravAcc(M1, P1, P2)
   T = P2-P1;
   T = -T*(6.67408*1e-11)*M1/(norm(T)^3);
end 

[Image Analyst]

What does this show in the command window:

whos S
whos AU
whos C1
whos C2
whos R
whos DT
bi

Chances are your index is more than the length of them, like S does not have bi+2 elements in it.

Answer 1

[Vimal Rathod]

You could try debugging your code by putting a breakpoint on the line of error and check if the concerned variables are defined and have expected values.