Propagating for an entire year is extremely difficult to get right. There are ways to do it, but you should not attempt to rederive or reimplement them yourself. Before I describe them, however, I have a different suggestion: don't propagate at all.
If you are designing a constellation of satellites that have a mission to perform, then the operators of those satellites are going to make station-keeping maneuvers, spending fuel to maintain their assigned orbits. Trying to model those maneuvers in addition to the perturbing forces would be even harder, but you don't have to, because the whole point of making the maneuvers is to cancel out the effect of the forces.
For example, due to the lumpiness of the Earth, satellites in geostationary orbits gradually drift east or west in subpoint longitude, but nobody lets that add up for a whole year unless the satellite is dead. In order to stay in their assigned boxes, all active, so-called stationary satellites regularly (roughly weekly) execute short burns (just a few seconds at a time is plenty) in order to correct that drift and move from the "downhill" edge of their box back to the "uphill" edge, and then start sliding down again until it comes time for the next burn back up the potential slope.
If you're trying to calculate what this looks like for a year, then your most accurate answer will be obtained by figuring out how many days of drift you have between burns, and just replaying those few days of motion enough times to fill out a year. Similar tricks are useful for low-orbiting satellites: find how many days it takes to nearly repeat a ground track, and simulate only that many days. This may not be an integer, and is unlikely to evenly divide 365, so you will actually get a better answer for coverage if you simulate only one repeat period than exactly a year.
Doing that to high precision requires some propagation, but not for a year, which you should work hard to avoid. Under many circumstances, though, even that much isn't necessary. For basic constellation design, figuring out the inclinations and altitudes of the orbit planes and how many satellites to put in each is the main thing, and it can be done reasonably well just from nominal orbits, as in a Walker circular layout.
If, however, you are stuck simulating objects which aren't actively controlled (whether you are tracking debris, or have a customer who for some reason is willing to pay for east-west burns but not for north-south burns to control inclination), then you have to do what you actually asked about. In that case, the core idea is to avoid position and velocity. They change much too fast, so your integrator will need a very small step size. This will force it to take so many steps that it will both take far too long and also produce a very incorrect answer. You need to make your integrator's life easier by enabling it to take many fewer, and thus much larger, steps. To maintain accuracy, you need to switch to different coordinates, that change much more slowly.
This is pretty much what Kepler's orbital elements already are. The reason we always refer to Kepler is that even with every imaginable perturbation turned on, most orbits don't vary much from an ellipse, at least for a while. I apologize for the way this becomes a link-only answer, but I am not writing an entire textbook in orbital mechanics. Even if I were, I still wouldn't try to shoehorn it all into one post. :) At some point, you just have to hit the library and read the back issues. The way to proceed is to look at Lagrange's planetary equations showing how the inclination, eccentricity, and so on change over time, as described here and here and here.
Once you've gone that far, stop to wonder whether there might be some other set of coordinates, in which the variation from long-term periodic motion is even easier to analyze. The answer is yes, and there are actually quite a few different kinds, but they get progressively more difficult even to describe, much less understand, derive, or implement well. Now take a look at Delaunay variables, and Brouwer's and Kozai's mean elements, and try to follow the path that led from there to the SGP4 we all love to hate, through the Goddard Trajectory Determination System (here, here, and here) to the Draper Semi-Analytic Satellite Theory (here and here), including developments published just two weeks ago.
Finally, note that SST, which has been under development for decades and published in dozens of journal articles and doctoral theses, is now part of the open-source orbit toolset called OreKit. Download that, and have a go at wrapping a Java library to leverage their already tested and peer-reviewed code, rather than trying to type in literally hundreds of pages of equations and somehow finding enough time or assistants to test them at all adequately.
