A big problem we get around $(x,y,z)=(0.822,1.265,1.855)$.
The Buffalo Way helps:
Let $x=\min\{x,y,z\}$, $y=x+u$,$z=x+v$ and $x=t\sqrt{uv}$.
Hence, $\frac{13}{5}\prod\limits_{cyc}(8x^3+5y^3)\left(\sum\limits_{cyc}\frac{x^4}{8x^3+5y^3}-\frac{x+y+z}{13}\right)=$
$$=156(u^2-uv+v^2)x^8+6(65u^3+189u^2v-176uv^2+65v^3)x^7+$$
$$+2(377u^4+1206u^3v+585u^2v^2-1349uv^3+377v^4)x^6+$$
$$+3(247u^5+999u^4v+1168u^3v^2-472u^2v^3-726uv^4+247)x^5+$$
$$+3(117u^6+696u^5v+1479u^4v^2+182u^3v^3-686u^2v^4-163uv^5+117v^6)x^4+$$
$$+(65u^7+768u^6v+2808u^5v^2+2079u^4v^3-1286u^3v^4-585u^2v^5+181uv^6+65v^7)x^3+$$$$+3uv(40u^6+296u^5v+472u^4v^2-225u^2v^4+55uv^5+25v^6)x^2+ $$
$$+u^2v^2(120u^5+376u^4v+240u^3v^2-240u^2v^3-25uv^4+75v^5)x+$$
$$+5u^3v^3(8u^4+8u^3v-8uv^3+5v^4)\geq$$
$$\geq u^5v^5(156t^8+531t^7+2t^6-632t^5-152t^4+867t^3+834t^2+299t+40)\geq0$$
Done!
For example, we'll prove that $$6(65u^3+189u^2v-176uv^2+65v^3)\geq531\sqrt{u^3v^3},$$ which gives a coefficient $531$ before $t^7$ in the polynomial $156t^8+531t^7+2t^6-632t^5-152t^4+867t^3+834t^2+299t+40.$
Indeed, let $u=k^2v$, where $k>0$.
Thus, we need to prove that:
$$130k^6+378k^4-177k^3-352k^2+130\geq0$$
and by AM-GM we obtain: $$130k^6+378k^4-177k^3-352k^2+130=$$
$$=130\left(k^3+\frac{10}{13}k-1\right)^2+\frac{k}{13}(2314k^3+1079k^2-5576k+2600)\geq$$
$$\geq\frac{k}{13}\left(8\cdot\frac{1157}{4}k^3+5\cdot\frac{1079}{5}k^2+21\cdot\frac{2600}{21}-5576k\right)\geq$$
$$\geq\frac{k^2}{13}\left(34\sqrt[34]{\left(\frac{1157}{4}\right)^8\left(\frac{1079}{5}\right)^5\left(\frac{2600}{21}\right)^{21}}-5576\right)>0.$$
We'll prove that $$
2(377u^4+1206u^3v+585u^2v^2-1349uv^3+377v^4)\geq2u^2v^2,$$ for which it's enough to prove that:
$$377t^4+1206t^3+584t^2-1349t+377\geq0$$ or
$$t^4+\frac{1206}{377}t^3+\frac{584}{377}t^2-\frac{1349}{377}t+1\geq0$$ or
$$\left(t^2+\frac{603}{377}t-\frac{28}{29}\right)^2+\frac{131015t^2-69589t+9633}{142129}\geq0,$$ which is true because
$$69589^2-4\cdot131015\cdot9633<0.$$
Since $e^x\geq 1+x$ for all $x\in \mathbb{R}$, by letting $x=-\log u$ we have that for any $u >0$, $\log u\geq 1-\frac{1}{u}$, and so $u\log u>u-1$ for any $u>0$. Combining these two inequalities, it follows that $$\sqrt{x^x y^y}= \exp\left(\frac{x\log x}{2}+\frac{y\log y}{2}\right)\geq \exp\left(\frac{x+y}{2}-1\right)\geq \frac{x+y}{2}.$$ Hence $$\frac{x^x}{x+y}+\frac{y^y}{y+z}+\frac{z^z}{z+x}\geq \frac{1}{2} \left(\sqrt{\frac{x^x}{y^y}}+\sqrt{\frac{y^y}{z^z}}+\sqrt{\frac{z^z}{x^x}}\right),$$ and the result follows from the AM-GM inequality.
Remark: By combining $e^x\geq 1+x$ for $x\in\mathbb{R}$ and $u\log u\geq u-1$ for $u>0$ as above, we obtain the following "reverse exponential AM-GM inequality": $$\sqrt[n]{x_{1}^{x_{1}}x_{2}^{x_{2}}\cdots x_{n}^{x_{n}}}\geq \frac{x_1+x_2+\cdots+x_n}{n}.$$
Best Answer
The standard way of solving the problem on a conditional extremum is the method of Lagrange multipliers, which reduces it to a system of equations.
The greatest value of function $$f(x,y,z,\lambda) = x^2y+y^2z+z^2x+\lambda(x+y^2+z^3-1)$$ on the interval $$x,y,z\in[0,1]$$ is reached or at its edges, or in the inner stationary point.
$\color{brown}{\textbf{Inner stationary points.}}$
The inner stationary points has zero partial derivatives $$\begin{cases} f'_\lambda = x + y^2 + z^3 - 1 = 0\\ f'_x = z^2 + 2xy + \lambda = 0\\ f'_y = x^2 + 2yz + 2\lambda y = 0\\ f'_z = y^2 + 2zx + 3\lambda z^2 = 0. \end{cases}$$ After the excluding of parameter $\lambda$ get the system $$\begin{cases} x + y^2 + z^3 - 1 = 0\\ x^2 + 2yz = 2y(z^2 + 2xy)\\ y^2 + 2zx = 3z^2(z^2 + 2xy), \end{cases}$$ or $$\begin{cases} x + y^2 + z^3 - 1 = 0\\ (1-3y^2-z^3)^2-4y^4+2yz(1-z)=0\\ 2z(1-3yz)(1-y^2-z^3)+y^2-3z^4=0. \end{cases}$$
Using of Groebner basis allows to get the positive solutions $$ \genfrac{[}{.}{0}{0}{x\approx 0.16367,\quad y\approx 0.761982,\quad z\approx 0.634724,\quad f\approx 0.454882} {x\approx 0.558113,\quad y\approx 0.498412,\quad z\approx 0.578371,\quad f\approx 0.485622}. $$ Alternative way is shown below.
$\color{brown}{\textbf{The edges.}}$
The edges of the field are achieved when $x = 0$, $y = 0$ or $z = 0$.
Substitution $x = 0$ in the expressions for the partial derivatives leads to the system $$ \begin{cases} x=0\\ y^2+z^3 = 1\\ 2yz+2\lambda y = 0\\ y^2+3\lambda z^2 = 0 \end{cases} $$ with solutions $$ \genfrac{[}{.}{0}{} {x=0,\quad y=\sqrt{0.75},\quad z=\sqrt[3]{0.25}\approx 0.629991,\quad f\approx 0.47247} {x=0,\quad y=0,\quad z=1,\quad f=0}$$
Substitution $y = 0$ in the expressions for the partial derivative leads to the system $$ \begin{cases} y=0\\ x+z^3=1\\ z^2+\lambda = 0\\ 2zx+3\lambda z^2 = 0 \end{cases} $$ with solution $$x=0.6,\quad y=0,\quad z=\sqrt[3]{0.4}\approx 0.736806,\quad f\approx 0.32573.$$
Substitution $z = 0$ in the expressions for the partial derivative leads to the system $$ \begin{cases} z=0\\ x+y^2=1\\ 2xy+\lambda = 0\\ x^2+2\lambda y = 0 \end{cases} $$ with solutions $$ \genfrac{[}{.}{0}{} {x=0.8,\quad y=\sqrt{0.2}\approx 0.447214,\quad z=0,\quad f\approx 0.286217} {x=0,\quad y=0,\quad z=1,\quad f=0} $$
The values of function at the vertices of the area (unit parallelepiped) equal to zero.
So the greatest value approximately equals to $0.485622$. Given accuracy of calculations provides the inequality $$\boxed{x^2y+y^2z+z^2x<1/2.}$$
$\color{brown}{\textbf{System resolving, alternative way.}}$
Set condition allows us to reduce the problem to finding unconditional extremes of $$f(y,z)=(1-y^2-z^3)^2y+y^2z+z^2(1-y^2-z^3).$$ Necessary optimality conditions in the field have the form: $$\begin{cases} f'_y = (1-y^2-z^3)^2-4y^2(1-y^2-z^3)+2yz-2yz^2 = 0\\ f'_z = -6yz^2(1-y^2-z^3)+y^2+2z(1-y^2-z^3)-3z^4 = 0, \end{cases}$$ or $$\begin{cases} 5y^4+y^2(6z^3-6)+y(2z-2z^2)+(z^3-1)^2 = 0\\ 6y^3z^2+y^2(1-2z)+6y(z^3-1)z^2-5z^4+2z = 0. \end{cases}$$
If to consider the coefficient of the highest power of $y$ as denominator in equation, and the remaining coefficients - numerators, we get the above equations. Then we can subtract the second equation factor $y$ from the first and repeat subtraction with factor $1$, obtaining the system
$$\begin{cases} C_{2,2}(z)y^2 + C_{2,1}(z)y + C_{2,0}(z) = 0\\ 6y^3z^2+y^2(1-2z)+6y(z^3-1)z^2-5z^4+2z = 0, \end{cases}$$ where $$C_{2,2}(z) = 36z^7-36z^4+20z^2-20z+5,$$ $$C_{2,1}(z) = 18z^6+102z^5-30z^2,$$ $$C_{2,0}(z) = 36z^{10}-72z^7+50z^5+11z^4-20z^2+10z.$$
Thus, the order of the first equation in $y$ reduced from fourth to second. Likewise, lowering the order of the second equation by a first equation, we obtain
$$\begin{cases} C_{2,2}(z)y^2 + C_{2,1}(z)y + C_{2,0}(z) = 0\\ D_{2,2}(z)y^2 + D_{2,1}(z)y + D_{2,0}(z) = 0,\\ \end{cases}\qquad(1)$$ where $$D_{2,2}(z) = 180z^8+576z^7-72z^5-144z^4+40z^3-60z^2+30z-5,$$ $$D_{2,1}(z) = 180z^7-30z^6-30z^5-60z^3+30z^2,$$ $$D_{2,0}(z) = 180z^{11}-252z^8+100z^6-28z^5+25z^4-40z^3+40z^2-10z.$$
The system $(1)$ is a linear in the unknowns $y^2$ and $y$, so $$y^2=\dfrac{\Delta_2(z)}{\Delta_0(z)},\quad y=\dfrac{\Delta_1(z)}{\Delta_0(z)},\qquad(2)$$ where $$\Delta_0(z) = C_{2,2}(z)D_{2,1}(z) - C_{2,1}(z)D_{2,2}(z),$$ $$\Delta_2(z) = -C_{2,0}(z)D_{2,1}(z) + C_{2,1}(z)D_{2,0}(z),$$ $$\Delta_1(z) = -C_{2,2}(z)D_{2,0}(z) + C_{2,0}(z)D_{2,2}(z),$$ or $$\Delta_0(z) = 90z^{10}-828z^9-1662z^8-144z^7+396z^6+1028z^5-200z^4+180z^3-220z^2+10z,$$ $$\Delta_2(z) = -90z^{13}+540z^{12}+30z^{11}+234z^{10}-864z^9-290z^8+256z^7+74z^6+220z^5-160z^4+100z^3-50z^2,$$ $$\Delta_1(z) = 576z^{13}-1296z^{10}+90z^9+815z^8+444z^7-5z^6-580z^5+126z^4+10z^3+80z^2-30z-5.$$
From $(2)$ for $\Delta_0(z)\not=0$ should be $$\Delta_1^2(z) - \Delta_2(z)\Delta_0(z) = 0,$$ or $$331776z^{26}-1484892z^{23}-19440z^{22}+1233720z^{21}+3079404z^{20}+195732z^{19}-3189924z^{18}-3109428z^{17}+368233z^{16}+2734116z^{15}+1243978z^{14}-741000z^{13}-805907z^{12}+75696z^{11}+164040z^{10}+82560z^9-172194z^8+53440z^7-10290z^6+32440z^5-7460z^4-4400z^3+100z^2+300z+25 = 0.$$
The coefficients are particially calculated using the Mathcad package, and $\mathcal{polyroots}()$ function is also used, which calculates all the roots of the polynomial by the "accompanying matrix" method.
Calculating values $y$ with $(2)$ and $x,f$ by the formula $$x = 1-y^2-z^3,\quad f=xy^2+yz^2+zx^2$$ and checking them by substituting in the original system, we obtain the following stationary points with the positive coordinates: $$ (x,y,z,f)\in\left[\genfrac{}{}{0pt}{0} {(0.5581125,\ 0.4984120,\ 0.5783713,\ 0.4856221)} {(0.1636702,\ 0.7619816,\ 0.6347238,\ 0.4548812)} \right. $$