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.
$$
Only a partial answer.
Assuming $a\le b\le c$, then we have $0<a\le1$ and $1\le c<3$. Now we have two cases, $b\le1$ and $b>1$. The case $b\le1$ is easy to deal with.
Assuming $0<a\le b\le1\le c<3$, we have
\begin{align}
\left(\frac{a+1}{a+b} \right)\left(\frac{b+1}{b+c}\right)\left(\frac{c+1}{c+a} \right) \ge 1
&\iff (a+1)(b+1)(4-a-b)\ge(a+b)(3-a)(3-b)\\
&\iff {\left(2-a-b\right)} {\left(1-a\right)} {\left(1-b\right)}\ge0.
\end{align}
Now assuming $0<a\le1\le b \le c<3$. Below is not an answer but only an analysis.
We want to know why this case is difficult, and why $\frac25$ is important.
Using the series expansion of $\exp$, and let $A = \operatorname{diag}\left(\ln\left(\frac{a+1}{a+b} \right),\ln\left(\frac{b+1}{b+c}\right),\ln\left(\frac{c+1}{c+a} \right)\right)$, then we have
\begin{align}
LHS
&= \sum_{n=0}^\infty\left(\frac25\right)^n\frac{\operatorname{tr}A^n}{n!}\\
&= 3 + \frac25\operatorname{tr}A + \frac12\cdot \left(\frac25\right)^2\operatorname{tr}A^2+ \frac16\cdot \left(\frac25\right)^3\operatorname{tr}A^3 +R_3,
\end{align}
where
$R_3 = \sum_{n=4}^\infty\left(\frac25\right)^n\frac{\operatorname{tr}A^n}{n!} \ge 0$, since $e^x-(1+x+x^2/2+x^3/6)$ is positive for any $x\in\mathbb R$.
Then we can prove the inequality if we have
$$\frac25\operatorname{tr}A + \frac12\cdot \left(\frac25\right)^2\operatorname{tr}A^2+ \frac16\cdot \left(\frac25\right)^3\operatorname{tr}A^3\ge\!\!\!?\;0,$$
which can be simplified to
$$75\operatorname{tr}A + 15\operatorname{tr}A^2+ 2\operatorname{tr}A^3\ge\!\!\!?\;0.\tag{1}$$
Numerical results suggest that using the 3 first terms is enough to prove the inequality. Note that in the first case where $b\le1$, using the first term $\operatorname{tr}A$ is enough (and what we did in the first part is in fact proving $\operatorname{tr}A\ge0$), that's why that case is easy.
So,
Why the case $b\ge1$ is difficult?
Because we have 2 more terms, $\operatorname{tr}A^2$ and $\operatorname{tr}A^3$, to deal with.
Why $\frac25$ is important?
Because $\frac25$ gives the coefficients 75, 15, and 2, which makes $75\operatorname{tr}A + 15\operatorname{tr}A^2+ 2\operatorname{tr}A^3\ge0$.
Best Answer
$$\mathbf{\color{green}{New\ proof,\ version\ of\ 02.07.18}}$$
$$\mathbf{\color{brown}{Using\ inequality}}$$
Let us consider the inequality $$(1+t)^{-\alpha}\geq1-\alpha t\quad \text{ for } t>-1,\quad \alpha>0,\tag1$$ which works due to the Maclauin series of $$(1+t)^{-\alpha} = 1 - \alpha t+\dfrac12(-\alpha)(-\alpha-1)t^2+\dfrac16(-\alpha)(-\alpha-1)(\alpha-2)t^3+\dots$$
Using of the inequality $(1)$ in the form of $$\left(\dfrac{a+1}{a+b}\right)^a=\left(\dfrac{a+b}{a+1}\right)^{-a}= \left(1+\dfrac{b-1}{a+1}\right)^{-a}\ge1-\dfrac{a(b-1)}{1+a}=2-b+\dfrac{b-1}{1+a}\tag2$$ is correct, because for the issue conditions $$a>0,\ b>0,\ c>0,\quad a+b+c = 3\tag3$$ $$\dfrac{b-1}{a+1} > -1.$$
$$\mathbf{\color{brown}{Task\ transformation}}$$
Using the inequality $(2)$ with the constraints $(3),$ easily to get $$\left(\dfrac{a+1}{a+b}\right)^a + \left(\dfrac{b+1}{b+c}\right)^b + \left(\dfrac{c+1}{c+a}\right)^c \geq 3+ \dfrac{b-1}{1+a}+ \dfrac{c-1}{1+b}+ \dfrac{a-1}{1+c},$$ so it can be proved the stronger inequality than the issue one: $$\dfrac{b-1}{1+a}+ \dfrac{c-1}{1+b}+ \dfrac{a-1}{1+c}\ge0,$$ or $$(b^2-1)(1+c)+(c^2-1)(1+a)+(a^2-1)(1+b)\ge0.\tag4$$ The least value of $LSH(4)$ can be achieved in the stationary points or on the edges of the area $(3).$
$$\mathbf{\color{brown}{Stationary\ points}}$$ The stationary points of the function can be found using the Lagrange multiplyers method, as the stationary points of the function $$f(a,b,c,\lambda)=(b^2-1)(1+c)+(c^2-1)(1+a)+(a^2-1)(1+b)+\lambda(a+b+c-3),$$ by the solving of the systen $f'_a = f'_b = f'_c = f'_\lambda=0,$ or \begin{cases} c^2-1+2a(1+b)+\lambda=0\\ a^2-1+2b(1+c)+\lambda=0\\ b^2-1+2c(1+a)+\lambda=0\\ a+b+c=3.\tag5 \end{cases} Summation of $(5.1)-(5.3)$ gives $\lambda=-4,$ then \begin{cases} c^2+2a(1+b)=a^2+2b(1+c)=b^2+2c(1+a)=5\\ a+b+c=3, \end{cases} \begin{cases} c^2+2a(4-a-c)=a^2+2(3-a-c)(1+c)=5\\ a+b+c=3, \end{cases} \begin{cases} 3c^2=3a^2-10a+4c+6\\ 3a^2-10a+4c+6+6a(4-a-c)=15\\ a+b+c=3, \end{cases} \begin{cases} (4-6a)c=3a^2-14a+9\\ 3(3a^2-14a+9)^2=(4-6a)^2(3a^2-10a+6)+4(4-6a)(3a^2-14a+9)\\ a+b+c=3, \end{cases} \begin{cases} (a-1)(27a^3-81a^2+45a+1)=0\\ c=\dfrac{3a^2-14a+9}{4-6a}\\ a+b+c=3,\tag6 \end{cases} The root $a=1$ leads to the solution $$a=b=c=1,\quad f(a,b,c,\lambda)=0.\tag7$$ Taking in account the task symmetry by the variables $a,b,c,$ the roots of cubic part of $(6.1)$ are the values of $a,b,c.$ At the same time, production of this roots is negative due to the Vieta theorem. This means that the solution $(7)$ is the single one, which satisfies conditions $(3).$
$$\mathbf{\color{brown}{The\ edges}}$$ To analyze the edges of area let us consider the case $a\to 0.$ It leads to the inequality in one variable $$(b^2-1)(2-b)+(3-b)^2-b-2\ge0,$$ or $$3-b(b-2)(b-3)\ge0.\tag8$$ Taking in account that for $b\in[2,3]$ $LSH(8)\ge3,$ and for $b\in(0,2)$ $$LSH(8)=3-b(2-b)(3-b)\ge3-\dfrac{(5-b)^3}9\ge0,$$ the inequality $(8)$ holds in the edges.
$$\mathbf{\color{green}{Proved}}$$
$$\mathbf{\color{black}{Old\ proof}}$$
$$\mathbf{\color{brown}{Task\ transformation}}$$
First, we should consider the inequality within the area. Using evident inequality $$\dfrac1{1-t}\geq1+t\quad \text{ for } t\in\left(-1,1\right)\qquad(1)$$ and AM-GM for $a,b,c>0,$ one can get: $$\dfrac{a+1}{a+b}=\dfrac{a+1}{3-c}=\dfrac12\,\dfrac{a+1}{1-\dfrac{c-1}2}\geq\dfrac{a+1}2\left(1+\dfrac{c-1}2\right) = \dfrac{a+1}2\,\dfrac{c+1}2\geq\sqrt{ac},$$ $$\left(\dfrac{a+1}{a+b}\right)^a + \left(\dfrac{b+1}{b+c}\right)^b + \left(\dfrac{c+1}{c+a}\right)^c \geq (ac)^{a/2} + (ba)^{b/2} +(cb)^{c/2}\\ =\exp\left(\frac{a}2\,\ln(ac)\right) + \exp\left(\frac{b}2\,\ln(ba)\right) + \exp\left(\frac{c}2\,\ln(cb)\right)\\ \geq 3\exp\left(\frac{a}6\,\ln(ac) + \frac{b}6\,\ln(ba)+\frac{c}6\,\ln(cb)\exp\right)\\ =3\exp\left(\frac{a+b}6\ln{a}+\frac{b+c}6\,\ln{b}+\frac{c+a}6\,\ln{c}\right) \\ =3\exp\left(\frac{3-c}6\,\ln{a}+\frac{3-a}6\,\ln{b}+\frac{3-b}6\, \ln{c}\right),$$ and that gives the possibility to prove the inequality $$(3-c)\ln{a}+(3-a)\ln{b}+(3-b)\ln{c}\geq0\tag2$$ for $a,b,c>0,\ a+b+c=3.$
$$\mathbf{\color{brown}{Stationary\ points}}$$ To do this, it suffices to find the least value of the function $$f(a,b) = (a+b)\ln{a}+(3-a)\ln{b}+(3-b)\ln(3-a-b).\tag3$$ The stationary points of the function can be found by equating to zero the partial derivatives $f'_a$ and $f'_b,$ or $$\begin{cases} \ln{a}+\dfrac{a+b}a -\ln{b} - \dfrac{3-b}{3-a-b} = 0\\[4pt] \ln{a}+\dfrac{3-a}b - \ln(3-a-b) - \dfrac{3-b}{3-a-b} =0, \end{cases}$$ $$\begin{cases} \ln{\dfrac{a}b} + \dfrac{a}{b} - \dfrac{a}{3-a-b} = 0\\[4pt] \ln{\dfrac{a}{3-a-b}} + \dfrac{3-a-b}{b} - \dfrac{a}{3-a-b} = 0,\\ \end{cases}$$ or $$\begin{cases} \ln{xy} - x + xy =0\\[4pt] \ln{x} -x + y =0\\[4pt] x=\dfrac{a}{3-a-b}\\[4pt] y=\dfrac{3-a-b}{b}, \end{cases}$$ $$\ln{y}+y(x-1)=0,\quad y=x-\ln{x},$$ $$\ln(x-\ln{x})+(x-1)(x-\ln{x})=0,\tag4$$ with single solution $$x=y=a=b=1,$$ which corresponds to the maximum $f(a,b)$ $$f_m=0.$$
$$\mathbf{\color{brown}{The\ edges}}$$ To analyze the edges of area let us consider the case $a=0.$ It leads to the inequality in
one variable $$\left(\dfrac{b+1}3\right)^b+\left(\dfrac{4-b}{3-b}\right)^{3-b}\geq 2,\tag5$$ which is satisfied for $0<b<3.$
Thus, $$\boxed{\left(\dfrac{a+1}{a+b}\right)^a + \left(\dfrac{b+1}{b+c}\right)^b + \left(\dfrac{c+1}{c+a}\right)^c \geq 3}.$$