N-Body Solar System Simulator - Why are there innacuracies in the x plane but not the others?

I have written a basic N-Body simulator that simulates the motion of the planets in the solar system.

The system reads the positions and velocities of the planets at a given time, measures the time for one full orbit of each planet and then compares the accuracy of the simulation against a highly accurate model.

I have analysed the data in the x, y and z planes. Interestingly, the simulation is highly accurate in the y and z planes but about 10x less accurate in the x plane. See attached graph.

Any idea why this is? I am thinking it is perhaps because the angular momentum from planets spinning on their axis has some effect on their orbit in the x plane? Or perhaps the effect of dark matter?

Here is the bit of the code that does the 'simulating'.

for nk = 2:length(LY)

t = 0;
dt = LY(nk)/num_it; 
G = 6.67e-11;
t_finish = LY(nk); %One orbital period for selected planet                      
N_t = round(t_finish/dt);
N = length(p); %find number of objects in simulation

%view in center of mass frame
v = v - mean((m*[1 1 1]).*v)/mean(m);

%find accelleration
a = a_v(p, m, G, sft);

%find starting energy
[KE,PE] = E(v,p,G,m);

%Main Simulation
for n = 1:N_t
    
    %'kick drift kick' second order technique is used here. This teqnique is
    %used as it preserves energy fairly efficiently. Source: https://en.wikipedia.org/wiki/Leapfrog_integration
    
    %first 1/2 kick
    v = v + a*dt/2;
    
    %drift
    p = p + v*dt;
    
    %new accellerations
    a = a_v(p,m,G,sft);
    
    %second 1/2 kick
    v = v + a*dt/2;
    
    %new time
    t = t + dt;
end
psave(nk,:) = p(nk,:);
vsave(nk,:) = v(nk,:);
end