I want to show that the first variational equation of the system is
$$\dot{\phi_1} = 1 +\frac{2x_0}{1-x_0t}\phi_1(t),\qquad\phi_1(0)=0$$
And solve the IVP approximately up to Order $O(\varepsilon^2)$. For the second point I have the hint to use a particular solution of the form $x_p(t)=a(1-x_0t)$.
I already got stuck on the first variational equation, which is defined in my script (Ordinary Differential Equations and Dynamical Systems, Teschl, p.46) as $\dot{y}=A(t,x)y$, where $A(t,x):=\frac{\partial f}{\partial x}(t,\phi(t,t_0,x))$.
I already have found the solution $$x(t) = \sqrt{\varepsilon}\tan\left(\sqrt{\varepsilon}t + \arctan\left(\frac{x_0}{\sqrt{\varepsilon}}\right)\right)$$
I would have assumed (as the derivative of $f(x)=x^2+\varepsilon$ is $f(x)=2x$) that my first variational equation would just be
$$\dot{y} = 2\cdot\left(\sqrt{\varepsilon}\tan\left(\sqrt{\varepsilon}t + \arctan\left(\frac{x_0}{\sqrt{\varepsilon}}\right)\right)\right)y$$
I guess I am not understanding the definition of the first variational equation correctly. It would be great if someone could explain to me where I went wrong.
In speaking of perturbations, the basis solution is the solution of the unperturbed problem, that for $ε$. So that gives $$ x(t)=\frac{x_0}{1-x_0t}. $$
The next idea is that small perturbations in the equation result in small perturbations of solutions with identical initial conditions, at least for some reasonable time span. There are multiple ways to introduce a perturbation term, the easiest is $x_ε(t)=x(t)+εv(t)+O(ε^2)$. Alternatively one could use $x_ε(t)=x(t)(1+εv(t)+O(ε^2))$ or $x_ε(t)=x(t)/(1+εv(t)+O(ε^2))$ or similar.
In the linear case one gets $$ \dot x+ε\dot v=x^2+2εxv+ε+O(ε^2)\implies \dot v=2xv+1. $$ The last is a linear differential equation for $v$, it can be solved via integrating factor.
Generally, if $\dot x=f(x,ε)$, then the basis solution solves $\dot x=f(x,0)$ and the perturbation solution as above can be developed as $$ \dot x+ε\dot v=f(x+εv,ε)=f(x,0)+\partial_xf(x,0)ε+\partial_εf(x,0)ε+O(ε^2) $$ removing the basis solution gives for the linear perturbation term $$ \dot v=\partial_xf(x,0)v+\partial_εf(x,v). $$ One can give separate names to the (partial) Jacobi matrices.
Idle thought: The given equation can also be characterized as Riccati equation. As such, its solutions can be parametrized as $x=-u'/u$, $u(0)=1$, $u'(0)=-x_0$, resulting in $$ \ddot u+εu=0. $$ Solving this as perturbation problem for $u=u_0+εu_1+...$ gives $$ \ddot u_0=0\implies u_0(t)=1-x_0t,\\ \ddot u_1 = -u_0\implies u_1(t)=-\frac{t^2}2+\frac{x_0t^3}6 $$ so that after reconstruction $$ x_ε(t)=-\frac{-x_0-ε(t-\frac{x_0t^2}2)}{1-x_0t-ε(\frac{t^2}2-\frac{x_0t^3}6)} $$