jmc

First, we can factor out $xy,$ to get $$xy(x^3+x^2+x+1+y+y^2+y^3)=xy(x^3+y^3+x^2+y^2+x+y+1).$$We know $x+y=3.$ Let $p=xy.$ Then $$9=(x+y{)}^2=x^2+2xy+y^2=x^2+2xy+y^2,$$so $x^2+y^2=9-2p.$

Also, $$27=(x+y{)}^3=x^3+3x^{2}y+3x{y}^2+y^3,$$so $x^3+y^3=27-3xy(x+y)=27-9p.$

Thus, $$\begin{array}{rl}xy(x^3+y^3+x^2+y^2+x+y+1)& =p(27-9p+9-2p+3+1)\\ & =p(40-11p)\\ & =-11p^2+40p\\ & =-11{\left(p-\frac{20}{11}\right)}^2+\frac{400}{11}\\ & \le \frac{400}{11}.\end{array}$$Equality occurs when $xy=p=\frac{20}{11}.$ By Vieta's formulas, $x$ and $y$ are the roots of $$t^2-3t+\frac{20}{11}=0.$$The discriminant of this quadratic is positive, so equality is possible. Thus, the maximum value is $\boxed{\frac{400}{11}}.$