In the following code
varx=var('x y')
paramu=var('A α β u')
paramuc=tuple(list(paramu)+[u])
paramb=var('R p_x p_y')
varl=var('λ')
varg=tuple(list(varx)+[λ])
def Cobb_Douglas(x, y, A, α, β) :
return A*x^α*y^β
def expense(x,y,p_x,p_y) :
return p_x*x+ p_y*y
L(x,y,λ) = Cobb_Douglas(x,y,A, α, β) -λ*(expense(x,y,p_x,p_y)-R)
stationary_points = lambda f: solve([gi==0 for gi in f.gradient()],varg)
opt=flatten(stationary_points(L))
opt1= [x.full_simplify().factor().simplify_rational()for x in opt]
show(opt1)
it is evident that $x$, $y$ and $\lambda$ could be factorized further. But I have not founded the way.
