Peculiar Amplitude of Earth's Orbit in Z Axis I'm currently building an N-body simulator for my masters project and I'm using JPL's Horizons system to compare against.

I started with just a couple orbits and everything seemed to be working nicely (After fixing many mistakes implementing a Fourth Order Integrator) but after extending to 10+ orbits I notice that in the Z-axis the Earth's amplitude increases linearly, according to JPL, but my simulation has a constant amplitude. JPL's result is shown in red in the bottom left:

  

What is causing this? Currently I only have the Sun, Earth, Moon and Jupiter in the system, but I didn't even realise the amplitude increased like this? What am I missing or got wrong? 

(Also, the graph only shows 3 simulated orbits, and 80 JPL orbits)
I can reproduce your JPL Z coordinate plot if I use Earth's heliocentric position relative to the J2000 ecliptic.
Naturally the amplitude has a minimum around the year 2000:



Relative to the ecliptic of date, which accounts for precession, the heliocentric Z coordinate of the Earth-Moon barycenter is much smaller:



The IAU 2006 precession model has two components, detailed in Capitaine et al. 2003:


precession of the ecliptic, a 47"/century change in the plane of Earth's orbit
precession of the equator, a 1.4&deg;/century change in the plane and axis of Earth's rotation


Fifty years of ecliptic precession amount to 23.5" or 0.000114 radian, and the first plot is consistent with that.
The classical precession of the equinoxes (the intersection of the two planes) is a combination of both components.

Williams 1994 examines Earth's equatorial precession (dominated by the Sun and Moon) and finds that Venus's effect is larger than Jupiter's.
When looking at fluctuations in Earth's year length, I noticed a cycle similar to Venus's synodic period.
I expect that adding Venus to your system would reproduce some ecliptic precession.