Find maximum of $xy$ subject to $(x+1)^2 + y^2 = 4$ without using Calculus techniques

If you are deadset on a way without calculus, lets start with your idea $2$, i.e. how to maximize $$p(x) = x^2y^2 = x^2\left(4-(x+1)^2 \right)= x^2(1-x)(3+x)$$

Now it is evident that this is a downward facing quartic with roots (with multiplicity) $-3, 0, 0, 1$, so the function is positive when $x \in (-3, 0)$ or $(0, 1)$. It must have maxima in those intervals and one of them must be the global maximum we seek.

Case $1$: $x \in (-3, 0)$
$$p = (-x)(-x) (1-x)(3+x)$$ is a product of four positive terms, so by AM-GM, if the sum were a constant this would be maximized when all terms are equal. Unfortunately the sum is not a constant, and the equality condition $(-x) = (1-x)$ is absurd. Both these problems however can be addressed by appropriate scaling. So consider maximizing for some real $\alpha, \beta > 0$, $$\alpha \beta p = (-x)(-x)(\alpha - \alpha x)(3\beta + \beta x)$$ Now the sum is a constant if we have $2+\alpha = \beta$ and for equality we must have $-x = \alpha -\alpha x = 3\beta + \beta x \implies \alpha = \frac14(\sqrt{33}-3), \beta = \frac14(5+\sqrt{33}), x = -\frac14(3+\sqrt{33})$

As all the conditions for AM-GM get satisfied here, the maximum in $(-3,0)$ is when $x = -\frac14(3+\sqrt{33}) = x_1$, say.

Case $2$: $x \in (0, 1)$
Following the same logic, we now have $$\alpha \beta p = (x)(x)(\alpha-\alpha x)(3\beta+\beta x)$$ and the conditions give $2+\beta = \alpha, x = \alpha - \alpha x = 3\beta + \beta x \implies x_2 = \frac14(\sqrt{33}-3)$

All that is left is to see which among $x_1, x_2$ gives you the maximum, and it is easily seen $p(x_1) > p(x_2)$, so the maximum is $\sqrt{p(x_1)} = \frac14\sqrt{\frac32(69+11\sqrt{33})} \approx 3.52$.

Your constraint is a circle centered at $(-1,0)$ with radius $2$. And the level curves of $xy$ with positive level are hyperbolae in the first and third quadrant. Since the circle is centered left of the origin, visualizing this it should be clear that the maximum value of $xy$ happens when both $x$ and $y$ are negative.

Furthermore, the curve $xy=M$ (where $M$ is the maximal value subject to the constraint) must intersect the circle in a tangency instead of at crossings, or else a slightly larger $M$ would be obtainable. For these two curves to intersect with tangency, there must be a repeated root to the equation that you get by eliminating one variable.

That is, if you eliminate $y$ to get $$(x+1)^2+\left(\frac{M}{x}\right)^2=4$$ and then expand and clear $x$ from the denominator to get $$x^4+2x^3-3x^2+M^2=0$$ we need to have a polynomial equation with a doubled root. So for some $a,b,c$:

$$\begin{align} &x^4+2x^3-3x^2+M^2\\ &=(x-a)^2(x-b)(x-c)\\ &=(x^2-2ax+a^2)(x^2-(b+c)x+bc)\\ &=x^4-(2a+b+c)x^3+(a^2+2a(b+c)+bc)x^2-(a^2(b+c)+2abc)x+a^2bc\\ \end{align}$$

which implies the system $$\begin{cases} 2a+b+c&=-2\\ a^2+2a(b+c)+bc&=-3\\ a^2b+a^2c+2abc&=0 \end{cases}$$

In the last equation we know $a\neq0$ so we can divide by $a$. Then we can eliminate $c$ by solving for it in the first equation and substituting it into the other two. All this gives: $$\begin{cases} -4a-3a^2-2b-2ab-b^2&=-3\\ -2a-2a^2-4b-4ab-2b^2&=0 \end{cases}$$

Double the top and subtract from the second to get: $$6a+4a^2=6$$ from which we can solve for $a$ using the quadratic formula. Remember that $a$ is the doubled root, and so it is the same as the $x$-value where $xy$ is maximal. Also we know this number is negative, which helps to discern between the two quadratic solutions. We find that $x=a=-\frac14\left(3+\sqrt{33}\right)$ at the point where $xy$ is maximal.

Now we can find: $$\begin{align} y&=-\sqrt{4-(x+1)^2}\\ &=-\sqrt{4-\left(-\frac14\left(3+\sqrt{33}\right)+1\right)^2}\\ &=-\sqrt{4-\left(\frac14\left(1-\sqrt{33}\right)\right)^2}\\ &=-\sqrt{4-\frac1{16}\left(34-2\sqrt{33}\right)}\\ &=-\frac14\sqrt{30+2\sqrt{33}} \end{align}$$

And so the maximal value is $$M=xy=\frac{1}{16}\left(3+\sqrt{33}\right)\sqrt{30+2\sqrt{33}}$$

I submit visual evidence that this is the maximal value. A plot of the circle and the curve $xy=M$. A larger $M$, and the two curves would not intersect.

Here is perhaps a shorter way to solve it without calculus. At an extremum we have a turning point, or the tangent is horizontal - so we may translate the graph to get a polynomial multiple of $x^2$. Note that if $p(x) = x^2(4-(x+1)^2)$, then

$$p(x+a)+C = x^2Q(x) - 2a(2a^2+3a-3)x+(C - a^2(1-a)^2(3+a))$$

now if there is an extremum, for some $a, C$ the polynomial above must have no linear or constant term. The constant term can always be set to zero by appropriate choice of $C$, but we will also need the linear term to have zero coefficient. $a=0$ corresponds to the known double root / minimum at $x=0$, hence we need

$$2a^2+3a-3 = 0 \implies a = -\tfrac14(3 \pm \sqrt{33})$$

Now all that is left is to test which root gives the larger value.