Problem with 4-body Matlab code

Question

I'm trying to model the 4 body problem to see how Jupiter, Earth and Mercury orbit the Sun. I found a two body script and adapted it as accordingly to modify my problem, but for some reason the plotting's created from my code do not seem to look as I feel that they should (the trajectories are rather chaotic) and I'm not really sure why...?

I have posted the code below along with the function that ode45 refers to, if anyone could help me find where I'm going wrong I'd be very grateful! In the code I have assumed that the masses are constant (being that of the sun) to make things easier but actual masses can be easily added in.

Thanks!!

Here is the code and the function:

 clear all
    % Set time interval of interest
    tspan = [0 1];    

    G = 5.951965229*(10^25); %Gravity of sun in AU^3/day^2;
    M=[1 9.543521107*(10^-4) 3.002730125*(10^-6) 1.660197494*(10^-7)]';     %masses in solar mass units

    %Initial positions and velocities

    %Initial sun position and velocity
    x00 = 0; y00 = 0; z00 = 0;
    vx00=0; vy00=0; vz00=0; 

    %Initial jupiter position
    x10=-5.290051789760673E+00;
    y10=1.201533314111204E+00; 
    z10=1.133818761444658E-01; 

    %Initial jupiter velocity
    vx10=-1.763381428668202E-03; 
    vy10=-7.009259824375218E-03;  
    vz10=6.853726794143213E-05; 

    %Intial earth position
    x20=-9.084371517338289E-01;  
    y20=3.932198307824979E-01; 
    z20=-1.106691686797854E-05;

    %Initial earth velocity
    vx20=-7.115556717957046E-03; 
    vy20=-1.584684870915936E-02;  
    vz20=1.736390456131702E-08; 

%Inital mercury position
    x30=-7.178277516331730E-03; 
y30=-4.625958019128616E-01; 
z30=-3.714027009954578E-02; 

%Inital mercury velocity
vx30=2.248843251006022E-02;  
vy30=1.006286364256332E-03; 
vz30=-1.980940891253320E-03; 

    z00=[x00 y00 z00 x10 y10 z10 x20 y20 z20 x30 y30 z30 vx00 vy00 vz00 vx10 vy10 vz10 ...
        vx20 vy20 vz20 vx30 vy30 vz30]'; %Initial z for 3d

    [t z] = ode45('fourbodyrhs', tspan, z00); 


    %The output is a column vector t, with a bunch of times,
    % and a matrix z.  Each row of z corresponds to one
    % time.  z has as many rows as t has elements.  Each row gives
    % the values of all of the components of z at one time.

    %Unpack the z matrix into variables you can understand
    x0 = z(:,1); y0 = z(:,2); z0 = z(:,3); x1 = z(:,4); 
    y1 = z(:,5); z1 = z(:,6); x2 = z(:,7); y2 = z(:,8);
    z2 = z(:,9); x3 = z(:,10); y3 = z(:,11); z3 = z(:,12); %3d

    %plot the trajectories of the two masses
    plot3(x0, y0, z0, x1, y1, z1, x2, y2, z2, x3, y3, z3)
    title('Four Body Trajectories')
    xlabel('x positions')
    ylabel('y positions')
    zlabel('z positions')

    axis equal

and here is the function 'fourbodyrhs'

function zdot =fourbodyrhs(t,z)

% This example is "The four body problem for planets orbiting the sun."

% Set constants
G = 5.951965229*(10^25); %Gravity of sun in AU^3/day^2
M=[1 9.543521107*(10^-4) 3.002730125*(10^-6) 1.660197494*(10^-7)]'; %masses in solar mass units


%Unpack the column vector z into things one can understand:
%four 2-element column vectors.
r{1} = z(1:3);    % The components of the position of sun
r{2} = z(4:6);    % The components of the position of jupiter
r{3} = z(7:9);    % The components of the position of Earth
r{4} = z(10:12);  % The components of the position of mercury
v{1} = z(13:15);  % Same as above but for velocities 
v{2} = z(16:18);    
v{3} = z(19:21);
v{4} = z(22:24);

N=4;

for i=1:N

    rdot{i} = v{i};
    vdot{i} = 0;

    for j=1:i-1

        vdot{i}=vdot{i}+G*M(j)*(r{j}-r{i})./((r{j}-r{i}).'*(r{j}-r{i})).^1.5;
    end

    for j=i+1:N

        vdot{i}=vdot{i}+G*M(j)*(r{j}-r{i})./((r{j}-r{i}).'*(r{j}-r{i})).^1.5;
    end

end


%Pack the 10 two-element column vectors above into one 20 element
%column vector
zdot= [ rdot{1}; rdot{2}; rdot{3}; rdot{4}; vdot{1}; vdot{2}; vdot{3}; vdot{4}];