Consider the nonlinear map $F_i:\mathbb R^2 \to \mathbb R$
$F_i(x):=\varepsilon^2\langle x, A_i x\rangle +\varepsilon\langle b_i,x \rangle + x_i,$
where $A_i$ is some matrix and $b_i$ some vector
Can we uniquely solve the equation $(F_1(x),F_2(x))=z_0 \in \mathbb R^2$ for $\varepsilon>0$ small enough such that $x = \sum_{n=0}^{\infty} \varepsilon^n z_n$ where $z_n$ are some coefficients?
This is a particular instance of the analytic implicit function theorem; see e.g. this MO page and further references there.
You can solve for the desired power series by starting with $x = x_0 + O(\varepsilon)$ and using iteration on the right-hand side of this form of your equation: $$ x = (I+\varepsilon b)^{-1} (x_0 - \varepsilon^2 A(x,x)) =: C(x), $$ where $(I+\varepsilon b)x = (x_1,x_2) + \varepsilon (\langle b_1, x\rangle, \langle b_2, x\rangle)$ and $A(x,y) = (\langle x, A_1 x\rangle, \langle x, A_2 x\rangle)$.
The question of convergence can be approached via the method of majorants (aka Cauchy majorants). Namely, first pick some matrix and vector norms so that $|b x| \le \|b\| |x|$ and $|A(x,x)| \le \|A\| |x|$, and pick some upper bounds $\|b\|<\beta$, $\|A\| < \alpha$. Then, comparing the coefficients of powers of $\varepsilon$, notice that $$ y = C(x) \quad \implies \quad \left\| y \right\| \preceq \frac{|x_0| + \varepsilon^2 \alpha \|x\|^2}{1+\varepsilon \beta} =: \Gamma(\|x\|) , $$ where $\|y\| = \sum_{n=0}^\infty \varepsilon^n |y_n|$ for $y= \sum_{n=0}^\infty \varepsilon^n y_n$ and $\preceq$ means that the the right-hand side is larger than the left-hand side on coefficient-by-coefficient basis as power series in $\varepsilon$ (the rhs majorizes the lhs). Of course, the power series $x$ converges iff the power series $\|x\|$ converges. By similar reasoning, $\|x\| \preceq \|y\|$ implies $\Gamma(\|x\|) \preceq \Gamma(\|y\|)$. So, defining $x^{(0)} = x_0$ and $x^{(k+1)} = C(x^{(k)})$, if we can find a sequence of convergent positive coefficient power series $\xi^{(k)} = \sum_{n=0} \varepsilon^n \xi_n^{(k)}$ such that $\xi^{(k+1)} \preceq \Gamma(\xi^{(k)})$ and it has a convergent limit $\xi^{(k)} \to \xi$, then we will have found a uniform majorant for each step of our iteration, $\|x^{(k)}\| \preceq \xi^{(k)} \preceq \xi$, meaning that our iteration also has a convergent limit $x^{(k)} \to x$. The easiest way to find such majorants is just to take a fixedpoint $\xi = \Gamma(\xi)$ and set $\xi^{(k)} := \xi$. As a quadratic equation, $\xi = \Gamma(\xi)$ has two solutions, but generally only one of them is a stable fixed point, which is given by $$ \xi = \frac{|x_0|}{1+\varepsilon \beta} \frac{2}{1 + \sqrt{1-\frac{4\varepsilon^2 \alpha |x_0|}{(1+\varepsilon\beta)^2}}} = \frac{|x_0|}{1+\varepsilon \beta} \frac{1}{1-\varepsilon^2 \alpha |x_0| + O(\varepsilon^2)} . $$ So the radius of convergence of $\xi$, and hence a lower bound on the radius of convergence of the solution $x$, can be estimated to be $\min \{ 1/\beta, (1/\sqrt{\alpha |x_0|})_* \}$, where by $({\cdots})_*$ we refer to the closest pole to the origin of the second fraction in the formula displayed above.
