$\renewcommand{\Re}{\operatorname{Re}}$ $\renewcommand{\Im}{\operatorname{Im}}$ $\newcommand{\bR}{\mathbb{R}}$ $\newcommand{\bC}{\mathbb{C}}$ $\newcommand{\bZ}{\mathbb{Z}}$ $\newcommand{\const}{\operatorname{const}}$
The most famous example are planets rotating around star which mass is much larger than masses of planets. Planets are attracted to the star and to one another according to Newton's gravity law.
If masses of planets are $0$ then movement of each planet is described by a separate ODE \begin{gather} \ddot{\boldsymbol{r}}_j= - (\boldsymbol{r}_j -\boldsymbol{r}_0)|\boldsymbol{r}_j -\boldsymbol{r}_0|^{-3}\qquad j=1,\ldots,n,\label{eq-4.A.1}\\ \boldsymbol{r}_0=0\label{eq-4.A.2} \end{gather} where $\boldsymbol{r}_0$ is the position of the star and mass of the star is $1$. Here and below $\dot{x}=\frac{d x}{dt}$, $\ddot{x}=\frac{d^2 x}{dt^2}$.
These equation one can integrate and derive Kepler laws; see Kepler laws of planetary motion.
On the other hand, if we do not neglect masses of the planets we get a system (in the coordinate system where the center of masses is at $0$) \begin{multline} \ddot{\boldsymbol{r}}_j = - (\boldsymbol{r}_j -\boldsymbol{r}_0)|\boldsymbol{r}_j -\boldsymbol{r}_0|^{-3} - \sum_{1\le k\le n, k\ne j} \mu_k (\boldsymbol{r}_j -\boldsymbol{r}_k)|\boldsymbol{r}_j -\boldsymbol{r}_k)|^{-3}, \\\ j=1,\ldots,n, \label{eq-4.A.3} \end{multline} \begin{gather} \boldsymbol{r}_0 = -\sum _{1\le k\le n}\mu_k \boldsymbol{r}_k \label{eq-4.A.4} \end{gather} where $\mu_1,\ldots, \mu_n$ are masses of planets relative to the sun.
In the general case we have many body problem which is really difficult even as $n=3$. However if $\mu_k\ll 1$ for all $k=1,\ldots, n$ then we are in the framework of regular perturbation theory and solution can be written as \begin{gather} \boldsymbol{r}_j = \sum_{\boldsymbol{\alpha}=(\alpha_1,\ldots,\alpha_n)} \boldsymbol{R}(t)_{j, \boldsymbol{\alpha}} \mu_1^{\alpha_1}\cdots \mu_n^{\alpha_n}. \label{eq-4.A.5} \end{gather}
Remark 1.
Here we have several small parameters and therefore we have a multipower series.
We assume that all distances (that is $\boldsymbol{r}_j - \boldsymbol{r}_k|$ with $(j,k,=0,1,\ldots, n$, $j\ne k$) are $\asymp 1$. Then one can prove that this asymptotic soluton works as \begin{gather} \mu |t| \ll 1, \qquad \mu:=\max (\mu_1,\ldots,\mu_n). \label{eq-4.A.6} \end{gather}
For larger times it fails due to secular motion (change of parameters of orbits with a speed $\asymp \mu $ .
So, there are two motions: slow secular motion and fast (actually, normal speed) periodic notion along orbits.
For real Solar system (see below) (\ref{eq-4.A.6}) is fulfilled for few hundreds of years. For longer periods one needs to use Multiscale Analysis; see Chapter 8.
Remark 2. Recurrent linear ODEs for $\boldsymbol{R}_{j, \boldsymbol{\alpha}}$ do not seem to be easy to integrate. However there are effective methods of integration because the original non-linear ODEs are not arbitrary ODEs but coming from Lagrangian mechanics; these calculations could be done using Hamiltonian mechanics. See f.e. PDE-textbook, Chapter 10.
Discussion 1. What does it mean for the real, rather an abstract Solar system? See f.e. Nasa table:
System Sun, Earth, Moon is described by \begin{align} &\ddot{\boldsymbol{r}}_1= -M\boldsymbol{r}_1|\boldsymbol{r}_1|^{-3} + m_2(\boldsymbol{r}_2-\boldsymbol{r}_1)|\boldsymbol{r}_2-\boldsymbol{r}_1|^{-3}, \label{eq-4.A.7}\\ &\ddot{\boldsymbol{r}}_2= -M\boldsymbol{r}_2|\boldsymbol{r}_2|^{-3} - m_1(\boldsymbol{r}_2-\boldsymbol{r}_1)|\boldsymbol{r}_2-\boldsymbol{r}_1|^{-3} \label{eq-4.A.8} \end{align} where $M$, $m_1$, $n_2$ are the mass of the Sun, $0$, $\boldsymbol{r}_1$, $\boldsymbol{r}_2$ are positions of Sun, Earth, Moon and we ignore the movement of Sun which is reasonable as $M\gg\gg m_1 \gg m_2$.
Then for $\boldsymbol{r}_*=(m_1+m_2)^{-1} \bigl( m_1 \boldsymbol{r}_1+ m_2 \boldsymbol{r}_2\bigr)$ the position of the center of mass of Earth and Moon and for $ \boldsymbol{r}= \boldsymbol{r}_2- \boldsymbol{r}_1$ we have equations \begin{multline} \ddot{\boldsymbol{r}}_* = - M \boldsymbol{r}_*| \boldsymbol{r}_*|^{-3} \\ + \underbracket{\frac{M m_1 \boldsymbol{r}_1 } {m_1+m_2} \bigl(|\boldsymbol{r}_1|^{-3}-|\boldsymbol{r}_*|^{-3}\bigr) +\frac{M m_2 \boldsymbol{r}_2 } {m_1+m_2} \bigl(|\boldsymbol{r}_2|^{-3}-|\boldsymbol{r}_*|^{-3}\bigr)} \label{eq-4.A.9} \end{multline} and \begin{gather} \ddot{\boldsymbol{r}}=-(m_1+m_2)\boldsymbol{r}|\boldsymbol{r}|^{-3} +\underbracket{M \bigl( \boldsymbol{r}_1 |\boldsymbol{r}_1|^{-3} -\boldsymbol{r}_2 |\boldsymbol{r}_2|^{-3}\bigr)} \label{eq-4.A.10} \end{gather} with \begin{gather} \boldsymbol{r}_1 = \boldsymbol{r}_* -\frac{m_2\boldsymbol{r}}{m_1+m_2},\qquad \boldsymbol{r}_2 = \boldsymbol{r}_* +\frac{m_1\boldsymbol{r}}{m_1+m_2}. \label{eq-4.A.11} \end{gather} Assuming that $|\boldsymbol{r}_*|\gg |\boldsymbol{r}|$ we see that \begin{gather*} |\boldsymbol{r}_j|^{-3} - |\boldsymbol{r}_*| ^{-3} = -3(\boldsymbol{r}_j - \boldsymbol{r}_*)\cdot \boldsymbol{r}_* |\boldsymbol{r}_*|^{-5} + O(|\boldsymbol{r}|^2 |\boldsymbol{r}_*|^{-5}). \end{gather*} Then the selected term in (\ref{eq-4.A.9}) is $O\bigl(M|\boldsymbol{r}|^2 |\boldsymbol{r}_*|^{-4}\bigr)$ while the selected term in (\ref{eq-4.A.10}) is $O\bigl(M|\boldsymbol{r}| |\boldsymbol{r}_*|^{-3}\bigr)$
Considering selected terms in (\ref{eq-4.A.9}) and (\ref{eq-4.A.10}) as perturbations we can employ usual series.
Discussion 2. What does it mean for real Sun, Earth, Moon system?