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.

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.
