The following response would be significantly shorter if you were allowed to use the Complex Analysis Theorem that
if $~\lim_{z \to z_0} a(z) = L \neq 0$
then $\lim_{z \to z_0} \frac{1}{a(z)} = \frac{1}{L}.$
I am assuming that the above theorem is out of bounds. Further, attempting to simply derive the above theorem (from scratch), would first require that some Real Analysis theorems be derived.
So, I will go for Kludginess.
This problem requires a backward and offbeat approach. A straightforward conversion of $z$ to $x + iy$ will be very messy. Fortunately, there is a middle ground.
Preliminary Concepts
(1)
Given $~r \in \Bbb{R}, ~s \in \Bbb{R^+},$
the constraint that $|r| < s$
is equivalent to the constraint that
$-s < r < s.$
(2)
Suppose that $u,v \in \Bbb{C}$
and that you want to establish that
$|u - v| < \epsilon$.
This may be established by showing that
$|\text{Re}(u - v)| < \frac{\epsilon}{2}~~~$ and
$~~~|\text{Im}(u - v)| < \frac{\epsilon}{2}.$
This is based on the idea that
$|\sqrt{\left[\frac{1}{2}\right]^2 + \left[\frac{1}{2}\right]^2} < 1.$
(3)
Suppose that $u,v \in \Bbb{C}$
and you are given that $0 < |u - v| < \delta.$
This implies that
$|\text{Re}(u - v)| < \delta~~~$ and
$~~~|\text{Im}(u - v)| < \delta.$
The justification for this is very similar to that applied in (2) above.
(4)
In this problem, rather than identifying a relationship between $\delta$ and $\epsilon$ such that $\delta = g(\epsilon)$, I will make things easy on myself by also requiring that $\delta \leq \frac{1}{10}$.
Thus, the relationship will look like $\delta = \min\left[\frac{1}{10}, g(\epsilon)\right].$
(5)
As suggested at the start of my answer, I am going to start with the requirement that
$\displaystyle \left| ~\frac{1}{z^2 + 1} - \frac{1}{2i + 1} ~\right| < \epsilon.$
Then, I will see where this leads in terms of the necessary constraint on $\delta$, which will always be required to be $\leq \frac{1}{10}.$
Assume that $z = (x + iy)$.
Let $w = (z^2 + 1)$ and let $\overline{w}$ represent the complex conjugate of $w$.
Let $f(z) = \displaystyle \left( ~\frac{\overline{w}}{|w|^2} - \frac{1 - 2i}{5} ~\right).$
It is desired that $|f(z)| < \epsilon.$
Based on (2) above, I will achieve this, by requiring that
$| ~\text{Re}[f(z)] ~| < \frac{\epsilon}{2}~~~$ and
$~~~| ~\text{Im}[f(z)] ~| < \frac{\epsilon}{2}.$
Using (1) above, this will be done by requiring that
$$-\frac{\epsilon}{2} < \text{Re}[f(z)] < \frac{\epsilon}{2} ~~~\text{and} ~~~ -\frac{\epsilon}{2} < \text{Im}[f(z)] < \frac{\epsilon}{2}. \tag{A} $$
The challenge is to specify $\delta$ so that
$$ 0 < \left| ~(x + iy) - (1 + i) ~\right| < \delta \leq \frac{1}{10} \tag{B} $$
implies that the constraints in (A) above are satisfied.
$w = (z^2 + 1) = (x^2 - y^2 + 1) + i(2xy).$
This implies that
$\overline{w} = (x^2 - y^2 + 1) - i(2xy)~~~$
and $~~~|w|^2 = (x^2 - y^2 + 1)^2 + (4x^2y^2).$
Using (1) and (3) above, (B) above implies that
$$1 - \delta < x < 1 + \delta ~~~\text{and}
~~~ 1 - \delta < y < 1 + \delta. \tag{C} $$
Since $~\delta \leq \frac{1}{10},~$ you have that
$(1 - \delta) > 0~~~$ and $~~~\delta^2 < \delta.$
Therefore,
$$ 1 - 2\delta < x^2 < 1 + 3\delta ~~~\text{and}
~~~ 1 - 2\delta < y^2 < 1 + 3\delta. \tag{D} $$
Consequently lower and upper bounds for $(x^2 - y^2 + 1)$ are
$$(1 - 5\delta) = (1 - 2\delta) - (1 + 3\delta) + 1 \leq (x^2 - y^2 + 1) \\ \leq [\text{similarly}] ~(1 + 5\delta). \tag{E} $$
Similarly, lower and upper bounds for $(-2xy)$ are
$$(-2 - 6\delta) < (-2)(1 + \delta)^2 < (-2xy) \\ <
(-2)(1 - \delta)^2 < (-2 + 4\delta). \tag{F} $$
Using (E) and (F), you have that
$$(1 - 5\delta) < \text{Re}(\overline{w}) < (1 + 5\delta) ~~~\text{and} \\ (-2 - 6\delta) < \text{Im}(\overline{w}) < (-2 + 4\delta). \tag{G}$$
Using (G) above, lower and upper bounds for $|w|^2$ are
$$(5 - 26\delta) < (1 - 5\delta)^2 + (-2 + 4\delta)^2 \leq |w|^2 \\ \leq (1 + 5\delta)^2 + (-2 - 6\delta)^2 < (5 + 95\delta). \tag{H}$$
Now, using (G), (H) above, lower and upper bounds for $~$ Re$\left(\frac{\overline{w}}{|w|^2}\right)~$
and $~$ Im$\left(\frac{\overline{w}}{|w|^2}\right)~$
may be identified.
A lower bound for $~$ Re$\left(\frac{\overline{w}}{|w|^2}\right)~$ is
$$\frac{1 - 5\delta}{5 + 95\delta} = \frac{1 + 19\delta}{5 + 95\delta} - \frac{24\delta}{5 + 95\delta} > \frac{1}{5} - \frac{24}{5}\delta > \frac{1}{5} - 5\delta.\tag{I} $$
An upper bound for $~$ Re$\left(\frac{\overline{w}}{|w|^2}\right)~$ is
$$\frac{1 + 5\delta}{5 - 26\delta} = \frac{1 - \frac{26}{5}\delta}{5 - 26\delta} + \frac{\frac{51}{5}\delta}{5 - 26\delta} < \frac{1}{5} + \frac{51}{10}\delta < \frac{1}{5} + 6\delta.\tag{J} $$
A lower bound for $~$ Im$\left(\frac{\overline{w}}{|w|^2}\right)~$ is
$$\frac{-2 - 6\delta}{5 - 26\delta} = \frac{-2 + \frac{52}{5}\delta}{5 - 26\delta} - \frac{\frac{82}{5}\delta}{5 - 26\delta} > -\frac{2}{5} - \frac{82}{10}\delta > -\frac{2}{5} - 9\delta.\tag{K} $$
An upper bound for $~$ Im$\left(\frac{\overline{w}}{|w|^2}\right)~$ is
$$\frac{-2 + 4\delta}{5 + 95\delta} = \frac{-2 - 38\delta}{5 + 95\delta} + \frac{42\delta}{5 + 95\delta} < -\frac{2}{5} + \frac{42}{5}\delta < -\frac{2}{5} + 9\delta.\tag{L} $$
Using (I) and (J)
$$-5\delta = \frac{1}{5} - 5\delta - \frac{1}{5} < \text{Re}[f(z)] < \frac{1}{5} + 6\delta - \frac{1}{5} = 6\delta.$$
Using (K) and (L)
$$-9\delta = -\frac{2}{5} - 9\delta + \frac{2}{5} < \text{Im}[f(z)] < -\frac{2}{5} + 9\delta + \frac{2}{5} = 9\delta.$$
$\delta$ must be chosen so that
$$-\frac{\epsilon}{2} < -5\delta, ~6\delta < \frac{\epsilon}{2}, ~-\frac{\epsilon}{2} < -9\delta, ~9\delta < \frac{\epsilon}{2}.$$
Coupled with the constraint that $\delta \leq \frac{1}{10}$
a workable final specification is
$$\delta = \min\left(\frac{1}{10}, \frac{\epsilon}{20}\right).$$