I am not sure why you have written the line
PlanetPositions = zeros(6,2);
so ode45 is not able to do anything. Also, it seems that you are trying to solve two independent systems of equations. Therefore, I suggest two calls to ode45. Run this code
%% Inputs
stepsize=1; %step size in days, should not exceed 1
finaltime=100; %days from t0
%% System Conditions
global sunm G
sunm=1.9891e30; %kg
G=1.4879*10^(-34); %AU^3/(kg*Day^2)
%% Sun Initial Conditions
x=0;
y=0;
z=0;
Vx=0;
Vy=0;
Vz=0;
SunInitial=[x,y,z,Vx,Vy,Vz]; %summary of initial conditions of the sun
%Note: not used in rest of code.
%% Earth Initial Conditions
x=-7.829690890802676E-01; %AU
y=6.009223815219861E-01; %AU
z=-1.377986550047871E-05; %AU
Vx=-1.075646922487205E-02;%AU/day
Vy=-1.370584048238274E-02;%AU/day
Vz=3.974096024543801E-08; %AU/day
EarthInitial=[x,y,z,Vx,Vy,Vz]; %summary of initial conditions of the earth
%% New planet
Mx=-2.829690890802676E-01; %AU
My=8.009223815219861E-01; %AU
Mz=-4.377986550047871E-05; %AU
MVx=-0.075646922487205E-02;%AU/day
MVy=-2.370584048238274E-02;%AU/day
MVz=5.974096024543801E-08; %AU/day
MercuryInitial=[Mx,My,Mz,MVx,MVy,MVz];
%% Call the Function
%outpu=call the function, time to run, initial conditions
[attime_earth,finalposition_earth]=ode45(@Propagator,0:stepsize:finaltime,EarthInitial);
[attime_mercu,finalposition_mercu]=ode45(@Propagator,0:stepsize:finaltime,MercuryInitial);
%% Plot the Result
hold on
plot3(finalposition_earth(:,4),finalposition_earth(:,5),finalposition_earth(:,6),'g'); %plot of position of earth in green
plot3(finalposition_mercu(:,4),finalposition_mercu(:,5),finalposition_mercu(:,6),'g--');
xlabel('Distance (AU)');ylabel('Distance (AU)');zlabel('Distance (AU)'); %label the axis
legend('Earth'); %tell me what the green line is
axis equal
hold off
%fix the axis so it has the same scale in x, y and z
function out = Propagator(t, PlanetPositions)
global sunm G
%state space representation idot=A1*PlanetPositions+A2(PlanetPositions)
%where idot is the derivative of the PlanetPositions
%this represents the linear part in A1 added to the non-linear part in A2
% |Vx| |0 0 0 1 0 0| |x | | 0 |
% |Vy| |0 0 0 0 1 0| |y | | 0 |
% |Vz| = |0 0 0 0 0 1| * |z | + | 0 |
% |Ax| |0 0 0 0 0 0| |Vx| | -G*sunm/R^2*|Rx| |
% |Ay| |0 0 0 0 0 0| |Vy| | -G*sunm/R^2*|Ry| |
% |Az| |0 0 0 0 0 0| |Vz| | -G*sunm/R^2*|Rz| |
% idot = A1 * PlanetPositions + A2(PlanetPositions)
% where |Rx| is the unit vector of R in the x direction
% where R=sqrt(x^2+y^2+z^2);
% PlanetPositions= zeros(6,2);
dx=zeros(6,1); %initialize the xdot matrix
A1=zeros(6,6);
A2=zeros(6,1);
A1(1,4)=1;
A1(2,5)=1;
A1(3,6)=1;
%% Earth
x=PlanetPositions(1); %x position
y=PlanetPositions(2); %y position
z=PlanetPositions(3); %z position
R = norm([x,y,z]); %find the radius from the earth to the sun
unitvector=[x,y,z]/R;
% do these calculations now so we only have to calculate it once
A2(4)=-G*sunm/R^2*unitvector(1); %acceleration in the x
A2(5)=-G*sunm/R^2*unitvector(2); %acceleration in the y
A2(6)=-G*sunm/R^2*unitvector(3); %acceleration in the z
TotAccelEarth=A1*PlanetPositions+A2; %sum up the linear + nonlinear accelerations
out=TotAccelEarth(:); %output the accelerations
end
