If $|b|$ and $|c|$ are to be as large as possible, it looks likely that $|a|$ will have to be small. So a first attempt would be to see what happens if $a=0$. In that case, the function becomes $|bz+c|$. That function takes the unit disc to a disc of radius $|b|$ with its centre at a distance $|c|$ from the origin. The condition for that image to lie in the unit disc is $|b| + |c| \leqslant 1$. The maximum value of $|bc|$ is then $\frac14$, occurring when $|b|$ and $|c|$ are both $\frac12$.
So my first guess was that the answer to the problem might be $\frac14$. But I then did some numerical experiments with graphs, and I found that if $a = 0.16$, $b = 0.56$ and $c = -0.56$, then the function $0.16z^2 + 0.56z - 0.56$ takes the unit disc to the green region in the diagram below, which clearly lies inside the (red) unit circle. But it has the property that $|bc| = 0.56^2 = 0.3136$, which is quite a bit larger than $\frac14$.
That is as far as I can go with this problem. I don't even have a guess as to what the maximum value of $|bc|$ might be.
This is a solution for the case where the parameters $a,b,c$ are all real. I suspect that the answer may be the same if they are allowed to be complex, but the proof would have to be quite different.
For the condition $|az^2 + bz + c| \leqslant 1$ to hold whenever $|z|\leqslant1$, it is sufficient (by the maximum modulus principle) to check that it holds when $|z|=1$. So we want to ensure that $|ae^{2i\theta} + be^{i\theta} + c| \leq 1$ (for all $\theta$). But if $a,b,c$ are real then $$\begin{aligned}|ae^{2i\theta} + be^{i\theta} + c|^2 &= (ae^{2i\theta} + be^{i\theta} + c)(ae^{-2i\theta} + be^{-i\theta} + c) \\ &= 2ac\cos(2\theta) + 2b(a+c)\cos\theta + a^2 + b^2 + c^2 \\ &= 2ac(2\cos^2\theta - 1) + 2b(a+c)\cos\theta + a^2 + b^2 + c^2 \\ &= 4acx^2 + 2b(a+c)x + (a-c)^2 + b^2, \text{ where }x = \cos\theta. \end{aligned}$$ So we want $f(x) = 4acx^2 + 2b(a+c)x + (a-c)^2 + b^2$ to have maximum value $1$ on the interval $-1\leqslant x \leqslant 1$. Differentiate to get $f'(x) = 8acx + 2b(a+c)$. That is zero when $x = -\dfrac{b(a+c)}{4ac}$, and at that point $$f(x) = \frac{b^2(a+c)^2}{4ac} - \frac{2b^2(a+c)^2}{4ac} + (a-c)^2 + b^2.$$ Put that equal to $1$ and solve for $b^2$, to get $$b^2 = 4ac\left(1 - \frac1{(a-c)^2}\right).$$ But for fixed $c$ we want to choose $a$ so as to maximise $b$. So differentiate $b^2$ with respect to $a$, getting $$4c\left(1 - \frac1{(a-c)^2}\right) + \frac{8ac}{(a-c)^3}.$$ Put that equal to zero to get $(a-c)^3 = -(a+c).$
Now let $t = a-c$. Then $t^3 = -a-c$, so that $c = \frac12(t^3 + t)$. Also, $t^6 - t^2 = (a+c)^2 - (a-c)^2 = 4ac$. Therefore $$b^2c^2 = (t^6 - t^2)\left(1 - \frac1{t^2}\right)\cdot\frac14(t^3+t)^2 = \frac14t^2(t^2-1)^2(t^2+1)^3.$$ On the interval $-1\leqslant t\leqslant 1$, that function has a maximum value $\dfrac{27}{256}$ (attained when $t =\pm \dfrac1{\sqrt2}$). So the maximum value of $|bc|$ is $\dfrac{3\sqrt3}{16} \approx 0.3248.$
From that, it is easy to see that this maximum value occurs when $a = \dfrac{\sqrt2}8 \approx 0.177$, $b = \dfrac{\sqrt6}4 \approx 0.612$ and $c =
-\dfrac{3\sqrt2}8 \approx -0.53$. Using those numbers, the diagram from the previous comment looks like this, with the green image of the unit disc under the function $|az^2 + bz + c|$ fitting snugly inside the unit circle.
\begin{tikzpicture} [scale=6]\draw [help lines, ->] (-1.1,0) -- (1.1,0) ; \draw [help lines, ->] (0,-1.1) -- (0,1.1) ; \draw [red] circle (1) ; \fill [green, domain=0:6.285, samples=100] plot ({0.177*cos(2*\x r) + 0.612*cos(\x r) - 0.53}, {0.177*sin(2*\x r) + 0.612*sin(\x r)}); \end{tikzpicture}