Note that the factorization below
suggests a the simpler equation,
i.e. translating $(-2,-3)$ to the origin
$$ y'=(x+2)(y+3)-6 $$
and the initial condition $(0,-1)$ to $(2,2)$.
This, however, leaves $(y+3)'=y'$ unchanged:
$$
\eqalign{
y' &= xy - 6 \\
y' - xy &= - 6 \\
(y' - xy) \, e^{-x^2/2} &= - 6\,e^{-x^2/2} \\
(y \, e^{-x^2/2})' &= - 6\,e^{-x^2/2} \\
y \, e^{-x^2/2} &= - 6 \int e^{-x^2/2} dx \\
y \, e^{-x^2/2} &= c-6 \sqrt{\frac\pi2}\,
\text{erf}\frac{x}{\sqrt{2}}
\qquad \implies \quad c=y_0 \\
y &= e^{x^2/2} \left( y_0 - 3 \sqrt{2\pi}\,
\text{erf}\frac{x}{\sqrt{2}}\right) .
}
$$
Translating back, we get
$$
y = e^{(x+2)^2/2} \left( c - 3 \sqrt{2\pi}\,
\text{erf}\frac{x+2}{\sqrt{2}}\right)-3
$$
for the same constant $c$ which, however,
is no longer equal to $y_0$, but rather to
$$
\eqalign{
c &= (y_0+3) e^{-\frac12(x_0+2)^2}
+ 3\sqrt{2\pi}\,\text{erf}\frac{x_0+2}{\sqrt{2}}
\\&=2e^{-2}+3\sqrt{2\pi}\,\text{erf}\sqrt{2}
\approx 7.44839864640532
}
$$
which would give
$$
\eqalign{
y(0.1)
&= e^{1/200} \left( c - 3 \sqrt{2\pi}\,
\text{erf}\frac{2.1}{\sqrt{2}}\right)-3
\\&\approx -1.21143171766941
\,.
}
$$
If we had just let $u=x+2$ and $v=y+3$ to begin with,
we could have rewritten our differential equation as
$v'=uv-6$
which is easier to differentiate by hand anyway:
$$
\eqalign{
\\v^{(1)}&=\left(u\right)\cdot&v+\left(-6\right)
\\v^{(2)}&=\left(u^{2}+1\right)\cdot&v+\left(-6 \, u\right)
\\v^{(3)}&=\left(u^{3}+3\,u\right)\cdot&v+\left(-6\,u^{2}-12\right)
\\v^{(4)}&=\left(u^{4}+6\,u^{2}+3\right)\cdot&v+\left(-6\,u^{3}-30\,u\right)
\\v^{(5)}&=\left(u^{5}+10\,u^{3}+15\,u\right)\cdot&v
+\left(-6\,u^{4}-54\,u^{2}-48\right)
\\v^{(6)}&=\left(u^{6}+15\,u^{4}+45\,u^{2}+15\right)\cdot&v
+\left(-6\,u^{5}-84\,u^{3}-198\,u\right)
\\v^{(7)}&=\left(u^{7}+21\,u^{5}+105\,u^{3}+105\,u\right)\cdot&v
+\left(-6\,u^{6}-120\,u^{4}-522\,u^{2}-288\right)
\\v^{(8)}&=\left(u^{8}+28\,u^{6}+210\,u^{4}+420\,u^{2}+105\right)\cdot&v
+\left(-6\,u^{7}-162\,u^{5}-1110\,u^{3}-1674\,u\right)
\\v^{(9)}&=\left(u^{9}+36\,u^{7}+378\,u^{5}+1260\,u^{3}+945\,u\right)\cdot&v
+\left(-6\,u^{8}-210\,u^{6}-2070\,u^{4}-5850\,u^{2}-2304\right)
}
$$
especially if one notices the recursion
$$
\left.\matrix{
v^{(n)}=a_nv+b_n\\\\
a_1=u,~b_1=-6}\right\}
\quad\implies\quad
\left\{\matrix{
a_{n+1}=a_n'+u\,a_n\\\\
b_{n+1}=b_n'-6\,a_n\,,}\right.
$$
which I managed to calculate in sage:
var('a,b,u,v,z')
a(u) = u
b(u) = -6
for n in range(1,10):
# print n, a, b #Note: a & b are functions of u
# print n, a(u).simplify_full(), b(u).simplify_full()
print '\\\\v^{(%d)}&=\\left(%s\\right)\cdot v+\\left(%s\\right)' \
% (n, latex(a(u).simplify_full()), latex(b(u).simplify_full()))
# print 'v^{(%d)}=%d' % (n, 2*a(2)+b(2))
z(u) = a
a(u) = (diff(a,u) + u*z).simplify_full()
b(u) = (diff(b,u) - 6*z).simplify_full()
Using the initial condition $(u,v)=(2,2)$
corresponding to $(x,y)=(0,-1)$, we get
$$
v^{(1)}=-2\\
v^{(2)}=-2\\
v^{(3)}=-8\\
v^{(4)}=-22\\
v^{(5)}=-76\\
v^{(6)}=-262\\
v^{(7)}=-980\\
v^{(8)}=-3794\\
v^{(9)}=-15428
$$
so that, factoring out the $(-1)$,
I computed the following terms and partial sums:
t = [1, 2, 2, 8, 22, 76, 262, 980, 3794, 15428]
x = 0.1
y = 0
for n in range(10):
tn = t[n]*x^n/factorial(n)
y = y + tn
print '%-4d %-14g %-16.12g' % (n, tn, y)
0 1 1
1 0.2 1.2
2 0.01 1.21
3 0.00133333 1.21133333333
4 9.16667e-05 1.211425
5 6.33333e-06 1.21143133333
6 3.63889e-07 1.21143169722
7 1.94444e-08 1.21143171667
8 9.40972e-10 1.21143171761
9 4.25154e-11 1.21143171765
So the analytic and Taylor series solutions agree.
Best Answer
The relevant Taylor series for computing a square root can be given by Newton's binomial series. Let $g : [1,\infty) \to [1,\infty),g(x)=\lfloor x^{1/2} \rfloor$. (To deal with $\sqrt{x}$ for $x \in (0,1)$, write it as $((x^{-1})^{1/2})^{-1}$.)
Then
$$x^{1/2}=\left (g(x)^2+ \left (x-g(x)^2 \right ) \right )^{1/2}=g(x) \left ( 1 + \frac{x-g(x)^2}{g(x)^2} \right )^{1/2} \\ = g(x) \sum_{k=0}^\infty {1/2 \choose k} h(x)^k$$
where $h(x):=\frac{x-g(x)^2}{g(x)^2}$. This ${1/2 \choose k}$ is a generalized binomial coefficient, given by $\frac{\Gamma(3/2)}{k! \Gamma(3/2-k)}$.
Now $\Gamma(3/2)=\sqrt{\pi}/2$. For $k>1$, this $\Gamma(3/2-k)$ is a bit mysterious, but we may connect it back to positive arguments of the $\Gamma$ function using the reflection formula:
$$\Gamma(z)=\frac{\pi}{\Gamma(1-z) \sin(\pi z)}$$
so
$$\Gamma(3/2-k)=\frac{\pi}{\Gamma(k-1/2) \sin(\pi(3/2-k))}=(-1)^{k+1} \frac{\pi}{\Gamma(k-1/2)}.$$
Thus ${1/2 \choose k}=\begin{cases} 1 & k=0 \\ (-1)^{k+1} \frac{\Gamma(k-1/2)}{2 \sqrt{\pi} k!} & k>0\end{cases}.$
So the series as a whole is
$$x^{1/2}=g(x) \left ( 1 + \sum_{k=1}^\infty (-1)^{k+1} \frac{\Gamma(k-1/2)}{2 \sqrt{\pi} k!} h(x)^k \right ).$$
The main point of this little computation was for us to realize that this method can be pretty slow. The exact speed depends on how close $x$ was to a perfect square. If we ignore constant coefficients, asymptotically you get the worst case when $x=n^2+n$ for some integer $n$, and in this case $h(x)$ is around $1/n$, or around $1/\sqrt{x}$. Similarly, $\frac{\Gamma(k-1/2)}{k!}$ behaves like $k^{-3/2}$, so in this worst case scenario you are looking at terms behaving like $k^{-3/2} x^{-k/2}$. Since the terms alternate, this means the relative error is decaying more or less like $(x^{-1/2})^k$. This is linear convergence: you get about $\log_{10}(x)/2$ correct decimal digits in each step (noting that you started with one correct digit).
In contrast, Newton's method (which for square roots is also called the Babylonian method) converges quadratically (loosely speaking, the number of correct digits doubles each iteration). Moreover Newton's method doesn't require us to precompute this $g(x)$, nor do we need an estimate of $\pi$, nor do we need to evaluate all these $\Gamma(k-1/2)$'s.
Note that the Taylor series method can be improved a bit by alternatively looking at $g(x)=\lceil x^{1/2} \rceil$, which will result in better performance when $x$ is slightly below a perfect square. But this $n^2+n$ scenario will behave pretty much the same either way.