Sorry for the long code but look at the Hessian. It is obviously false since $H_{33} = 0$ not $-p_x$
%display latex var('A, x, y, l, alpha, beta, R, p_x, p_y'); U= Ax^(alpha)y^(beta); show(LatexExpr("U(x)="+latex(U))) D = p_xx + p_yy; show(LatexExpr("D="+latex(D))) L = U-l*(D-R) show(LatexExpr("L(x,y,l)="+latex(L))) L_x= L.diff(x) show(LatexExpr("L_x(x,y,l)="+latex(L_x))) L_y= L.diff(y) show(LatexExpr("L_y(x,y,l)="+latex(L_y))) L_l= L.diff(l) show(LatexExpr("L_l(x,y,l)="+latex(L_l))) z=solve([L_x==0, L_y==0, L_l==0,], x, y, l) show(z[0]) x1=z[0][0].right() show(LatexExpr("x^\star(p_x,p_y,R)="+latex(x1))) y1=z[0][1].right() show(LatexExpr("y^\star(p_x,p_y,R)="+latex(y1))) l1=z[0][2].right().canonicalize_radical() show(LatexExpr("l^\star(p_x,p_y,R)="+latex(l1))) U1=U.subs(x=x1,y=y1).canonicalize_radical() show(LatexExpr("U^\star(p_x,p_y,R)="+latex(U1))) L_xx= L_x.diff(x); L_xy= L_x.diff(y); L_xl= L_x.diff(l); L_yx= L_y.diff(x); L_yy= L_y.diff(y); L_yl= L_y.diff(l); L_lx= L_l.diff(x); L_ly= L_l.diff(y); L_ll= L_l.diff(l); H=matrix([[L_xl,L_xx, L_xy],[L_yl,L_yx, L_ly],[L_xl,L_lx, L_ly]]); H
Also if I typeset
show(LatexExpr("U^\star(p_x,p_y,R)="+latex(U1)))
I have only the Latex code, which is astonning since SageAMth is capable to print H