To answer your questions,
- Do the terms (J2, etc.) change with time?
Yes, they do. The Earth is still rebounding from the end of the last glaciation, and the Earth's rotation rate is decreasing due to transfer of angular momentum to the Moon's orbit. The end result is that $J2$ is decreasing by about 3 parts per million per century. You don't need to model that, however. The gravity models that GMAT uses are static models. You can enable tides, making the gravity model somewhat dynamic, but the dynamics are purely cyclical. There are no secular terms.
- I found other terms coefficients here. How to add them also to differential equation?
Find a software package that uses spherical harmonics to model gravitation. Do not roll your own. And if you do use such a package, make sure the coefficients you use are consistent with the model. The coefficients you found are denormalized. Most modern spherical harmonics models expect fully normalized coefficients.
A much trickier issue lies in the $\bar C_{20}$ term. Some of the Earth's equatorial bulge results from tidal interactions with the Moon and the Sun. The spherical harmonics gravitational coefficients are typically computed as if the Moon and Sun were not present. These tide-free models omit the contribution of the Moon and Sun to the equatorial bulge, and hence to $J_2$. (Technically, the frequency coefficients used to model the Earth tides have a zero frequency term that is nonzero.) The Earth spherical harmonics coefficients used by GMAT are tide-free. If you want to model how satellites precess due to the tidal bulge, but don't want to use a full blown Earth tide gravity model, it's better to add a tiny bit to the the $\bar C_{20}$ term. See Section 6, Gravitation of IERS Technical Note 36 for details.
- Why does the shown equation give 1 km error in 1 day in comparison with GMAT (configuration is below)?
The GMAT JGM-2 gravity model is in the file JGM2.cof of the GMAT source code tree. The first few non-comment lines in this file are
POTFIELD 70 70 1 3.98600441500000e+14 6.37813630000000e+06 1.00000000000000e+00
RECOEF 2 0 -4.84165390000000e-04
RECOEF 2 1 -1.86987640000000e-10 1.19528010000000e-09
RECOEF 2 2 2.43908370000000e-06-1.40010930000000e-06
You need to use compatible values to have your Julia integrator be consistent with GMAT's implementation of JGM-2. The first value, $3.986004415\times10^{14}$, is the TT-compatible value of the Earth's gravitational coefficient in $\text{m}^3/\text{s}^2$. You should be using this rather than the WGS-84 value of $3.986004418\times10^{14}$. The WGS-84 model is aimed at GPS; it is relativistically correct and hence uses a TDB-compatible value of $GM_\oplus$.
The second value, $6.3781363\times10^{6}$, is the equatorial radius of the Earth, in meters. You should be using this rather than the WGS-84 value if you want to try to match the results from GMAT.
The next three lines contain the normalized values of the cosine and sine coefficients of the JGM-2 gravitational potential model for the Earth. The first, -4.84165390000000e-04, is the tide-free value of the fully normalized $\bar C_{20}$ coefficient. The $\bar C_{20}$ is directly related to $J_2$ via $J_2 = -\sqrt{5}\bar C_{20}$.
From poking at the GMAT source code, there is a correction for the $\bar C_{20}$ term gravity coefficient if that term includes the permanent tide and if the user enables Earth tides. There does not appear to be a correction for a tide-free value that adds the permanent tide effect if the user disables Earth tides. The widely published value of $J_2$, 0.0010826359, includes the permanent tides. You probably should be using a value consistent with the tide-free value of $\bar C_{20}$, or $\sqrt{5}\,4.8416539\times10^{-4} \approx 0.00108262672$ instead of that widely published value.
The next two lines contain $\bar C_{21}$, $\bar S_{21}$, $\bar C_{22}$, and $\bar S_{22}$. Note that the $\bar C_{21}$ and $\bar S_{21}$ are nonzero. This means that the Earth's axis of rotation is not quite in line with the way we measure latitude and longitude. The Earth's rotation axis undergoes a small polar motion; you should turn this off in GMAT (if you can) to have a better chance of having your integrated state agreeing with the state computed by GMAT. If you can, you should also turn off Earth nutation and precession in GMAT.
Finally, you're using a different language, and possibly a different numerical integrator. You shouldn't be surprised on seeing differences with exactly the same code when one changes compiler optimization options or migrates to a different computer. You shouldn't be surprised at all when using a different language.