Wrong hessian

edit

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= A*x^(alpha)*y^(beta);
show(LatexExpr("U(x)="+latex(U)))
D = p_x*x + p_y*y;
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