I have a system of ODEs defined on $\mathbb{R}\times\mathbb{S}$, $$\begin{aligned}\dot{x}={ }&y\\\dot{y}={ }&-0.2y+\frac{300\cos(2)\sin(x)}{1.8(1.3+\cos(x-2))-2\sin(x)\sin(2)},\end{aligned}$$
This system has a saddle at $(0,0)$ and I want to compute the unstable manifold which is supposed to be pretty close to a homoclinic trajectory. So, I eliminate time, and solve the resulting DE \begin{equation}\frac{dy}{dx}=-0.2+\frac{1}{y}\frac{300\cos(2)\sin(x)}{1.8(1.3+\cos(x-2))-2\sin(x)\sin(2)}\tag{*}\end{equation} using power series $y(x)=a_1 x+a_2x^2+\dots+a_kx^k +\mathcal{O}(x^{k+1})$ for different $k$. That is, I plug my $y(x)$ into (*) and compute the coefficents $a_i$. I do so using symbolic mathematics, so there are shouldn't be any stupid mistakes (hopefully).
My problem: For all $k>3$, the solutions behave pretty well halfway to the point when it comes closer to the equilibrium. But then it diverge wildly from the expected trajectory. When I add more terms to the expansion (i.e., increase $k$), the divergence becomes even bigger. See the picture below.
Using numerical simulation I got a solution that leaves $(0,0)$ and makes an arc toward $(2\pi,0)$. In this sense, the best result seems to be given by either $k=3$ or $k=7$, but they are still very far from what I expect.
Following @LutzL's suggestion I tried to determine the radius of convergence of the power series expansion of the fraction in the second DE. The denominator does not have any real zeroes, so I computed a complex zero, but it does not restrict the radius of convergence: $x_0\approx -227.7655 + 3.15i$. However, I'm not sure if this pole is unique.
The problem was originally formulated for the general system \begin{align}\dot{x}={ }&y\\\dot{y}={ }&-ky+f(x),\end{align} with $2\pi$-periodic $f$, $f(2k\pi)=0$, $f'(0)>0$ so that the system has a saddle point at $(0,0)$. The solution curves as graphs of functions $y(x)$ also satisfy the DE $$ \frac{dy}{dx}=-k+\frac{f(x)}{y(x)} $$ The solutions with limit in the saddle point have initial condition $y(0)=0$, which is a singular point, and are determined by $y'(0)=r$ where $r$ is a solution of $r^2+kr=f'(0)$.
After multiplying with $y$ an alternative form of the equation is $$ \frac12\frac{d(y^2)}{dx} = -ky+f(x),~~ y(x)=\sqrt{2\int_0^x[-ky(s)+f(s)]ds} $$ which introduces a different kind of singularities (limiting the radius of converngence) at the roots of the integral on the right, which are branching points of the square root. At the same time they are also roots of $y$.
Due to the structure of that equation one has to expect trouble at the roots of $y$ and of course at the poles of $f$. Not all of them are seen in the numerical solution, but the power series is also valid as solution for the differential equation over the complex plane
The problem as specified now uses $k=0.2$ and $$ f(x)=\frac{300\cos(2)\sin(x)}{1.8(1.3+\cos(x-2))-2\sin(x)\sin(2)}. $$ This leads to an initial slope around $r=14$. To get rational coefficients to use in the CAS Magma which can deal in truncated power series with rational, but not with floating point coefficients, change $f$ to $$ f(x)=(r^2-k^2/4) \,\frac{\frac{33}{20} \, \sin(x)}{\frac{12}5-\frac34\cos(x)-\frac15\sin(x)} $$ where $r=14$ and $k=\frac15$. The modified coefficients are within $5\%$ of the original. Compute the power series for $y$ up to $O(x^{60})$ and with that compare the partial series to the numerical solution. Plot these in groups, here of 12 graphs. Alongside compute and plot the roots of the power series with a circle of radius $1.9$.