Maximizing $\frac{a^2+6b+1}{a^2+a}$, where $a=p+q+r=pqr$ and $ab=pq+qr+rp$ for positive reals $p$, $q$, $r$ 
Given $a$, $b$, $p$, $q$, $r \in\mathbb{R_{>0}}$ s.t.
  $$\begin {cases}\phantom{b}a=p+q+r=pqr \\ab =pq+qr+rp\end{cases} $$
  Find the maximum of $$\dfrac{a^2+6b+1}{a^2+a}$$


This question is terrifying that I even don't know how to start it. I've found different ways such as multiplying the numerator and denominator by $a$ and substitute, but it is not useful (maybe?). Anyway, please comment or answer if you solve it or having clues that may help solve this question.
 A: Hint.
Using the AM-GM we get
$$
a \ge 3\sqrt[3]{pqr}=3\sqrt[3]{a}\Rightarrow a\ge \sqrt{27}
$$
$$
ab\ge 3\sqrt[3]{a^2}\ge 3\sqrt[3]{27}
$$
now we can solve
$$
\max \frac{a^2+6b+1}{a^2+a}\ \ \ \text{s. t.}\ \ \ (a\ge \sqrt{27})\cap (ab \ge 9)
$$
and the solution is
$$
a = 3\sqrt 3,\ \ b = \sqrt 3
$$
with
$$
\frac{a^2+6b+1}{a^2+a} = 1+\frac{1}{3 \sqrt{3}}
$$
NOTE
Using the Lagrangian to obtain the stationary points
$$
L(a,b,\lambda,\mu,s_1,s_2) = \frac{a^2+6b+1}{a^2+a}+\lambda (a-\sqrt 27-s_1^2)+\mu(a b-9-s_2^2)
$$
the stationary points are the solutions for
$$
\nabla L = 0 \Rightarrow\left\{
\begin{array}{rcl}
 \frac{2 a}{a^2+a}-\frac{(2 a+1) \left(a^2+6 b+1\right)}{\left(a^2+a\right)^2}+\lambda +b \mu & = & 0 \\
 a \mu +\frac{6}{a^2+a} & = & 0\\
 -s_1^2+a-3 \sqrt{3} & = & 0 \\
 -s_2^2+a b-9 & = & 0 \\
 \lambda  s_1 & = & 0\\
 \mu  s_2 & = & 0 \\
\end{array}
\right.
$$
A: Let $p+q+r=3u$, $pq+pr+qr=3v^2$ and $pqr=w^3$.
Thus, $$w^3=3u,$$
$$a=3u,$$
$$b=\frac{v^2}{u}$$ and
$$\frac{a^2+6b+1}{a^2+a}=\frac{9u^2+\frac{6v^2}{u}+1}{9u^2+3u}=\frac{9u^3+6v^2\sqrt{\frac{w^3}{3u}}+\frac{w^3}{3}}{9u^3+3u^2\sqrt{\frac{w^3}{3u}}},$$ which is increasing function of $v^2$, which says that it's enough to find a maximal value of the last expression for equality case of two variables and since the last expression is symmetric of $p$, $q$ and $r$ and homogeneous  which degree $0$, it's enough to do it for $q=r=1.$
Can you end it now?
A: Here is an approach, which provides, so I think, a certain understanding of the "working domain", prior to the consideration of maximisation issue.
Your issue is to find values of $a$ and $b$ such that polynomial equation :
$$(x-p)(x-q)(x-r)=0 \ \ \iff$$
$$x^3-(p+q+r)x^2+(qr+pr+pq)x-pqr=0 \ \ \iff \tag{2}$$
$$x^3-ax^2+abx-a=0 \tag{1}$$


*

*(i) has three real roots ($p,q$ and $r$),

*(ii) that are all positive (implying that $a$ and $b$ are themselves positive).
Take a look at the following figure. Each little red circle represents a couple of (random) values $(a,b)$ ($a$ : abscissas, $b$ : ordinates) fulfilling condition (i). Among them, clearly, only the circles with a star also fulfill condition (ii), occupying a very tiny domain...

What is the limit of these circles, i.e., how has the green curve been obtained ? Its equation is $d=0$ where $d$ denotes the so-called discriminant $d$ of parametric equation (1), a particular case of a "resultant", concepts that you will see in University if you do a degree in mathematics. See remark 3 below. Here it is given under the form of a determinant : 
$$d=\begin{vmatrix}
1&-a&ab&-a&0\\
0&1&-a&ab&-a\\
3&-2a&ab&0&0\\
0&3&-2a&ab&0\\
0&0&3&-2a&ab\\
\end{vmatrix}=-a^2(a^2b^2-4a^2-4ab^3+18ab-27)\tag{3}$$
(maybe, you have recognized in the two first rows of (3) the entries of polynomial (1) and on the 3rd, 4th and 5th rows the entries of its derivative, another concept that hopefuly you haven't met yet. Please note the progressive shifting).
Now, the maximisation issue : one finds (I don't give a proof) that it is the point at the extreme left of the "good" spiky region that achieves the maximality, which is exactly point $(3\sqrt{3},\sqrt{3})$ found by  @Cesareo.
Remarks : 
1) Considering $d=0$ as a quadratic in variable $a$, with entries depending on parameter $b$, one can express $a$ as a function of $b$ :
$$a=\dfrac{2b^3-9b\pm\sqrt{\delta}}{b^2-4} \ \ \text{where} \ \ \delta=b^6-18b^4+108(b^2-1),$$
allowing to plot the green curve rather easily .
2) Formulas dealing with $p+q+r$, $pq+qr+rp$ and $pqr$ are called Vieta's formulas.
3) Ask your professor why a discriminant equal to $0$, already for a quadratic equation, expresses the fact that there are double roots.
A: Recently, I find a solution which is quite simple. Check it out!

Claim $1$: $a^2\ge3ab\rightarrow a\ge3b, b\le\dfrac{a}{3}$

Prove:
$$a^2-3ab=\left(x+y+z\right)^2-3\left(xy+yz+zx\right)=x^2+y^2+z^2-xy-yz-zx=\dfrac{\left(x-y\right)^2+\left(y-z\right)^2+\left(z-x\right)^2}{2}\ge0 \\ \therefore a^2\ge 3ab$$

Claim $2$: $a\ge3\sqrt{3}$

Prove:
By AM-GM Inequality,
$$\left(\dfrac{x+y+z}{3}\right)^3\ge xyz\\\left(\dfrac{a}{3}\right)^3\ge a\\\dfrac{a^3}{27}\ge a\\a^2\ge27\\\therefore a\ge3\sqrt{3}$$
Then comes the solution:
$\dfrac{a^2+6b+1}{a^2+a}\le\dfrac{a^2+2a+1}{a^2+a}=\dfrac{\left(a+1\right)^2}{a\left(a+1\right)}=\dfrac{a+1}{a}=1+\dfrac{1}{a}\le1+\dfrac{1}{3\sqrt{3}}=\boxed{1+\dfrac{\sqrt{3}}{9}}$
Equality holds when $x=y=z=b=\sqrt{3},a=3\sqrt{3}$
