from sympy import symbols, Symbol, Eq, Derivative, solve, Function, Rational, UnevaluatedExpr
from sympy.vector import Del

h = symbols('h') # labor in hours
s = symbols('s') # steel in tons
lamda = symbols('lamda') # steel in tons
R = Function('R')(h, s) #revenue
g = Function('g')(h, s) #constraint

Del_R = Symbol(r'\nabla R')
Del_g = Symbol(r'\nabla g')

R_eq = Eq(R, 100*h**Rational(2, 3)*s**Rational(1, 3))
R_eq

$\displaystyle R{\left(h,s \right)} = 100 h^{\frac{2}{3}} \sqrt[3]{s}$

B = 20000

g_eq = Eq(Eq(g, 20*h + 2000*s), B, evaluate = False)
g_eq

$\displaystyle g{\left(h,s \right)} = 20 h + 2000 s = 20000$

Lm_eq = Eq(Del_R, lamda*UnevaluatedExpr(Del_g))
Lm_eq

$\displaystyle \nabla R = \lambda \nabla g$

diff_R_eq_h = Eq(Derivative(R, h), Derivative(100*h**Rational(2, 3)*s**Rational(1, 3), h)).doit()
diff_R_eq_h

$\displaystyle \frac{\partial}{\partial h} R{\left(h,s \right)} = \frac{200 \sqrt[3]{s}}{3 \sqrt[3]{h}}$

diff_R_eq_s = Eq(Derivative(R, s), Derivative(100*h**Rational(2, 3)*s**Rational(1, 3), s)).doit()
diff_R_eq_s

$\displaystyle \frac{\partial}{\partial s} R{\left(h,s \right)} = \frac{100 h^{\frac{2}{3}}}{3 s^{\frac{2}{3}}}$

diff_g_eq_h = Eq(Derivative(g, h), Derivative(20*h + 2000*s, h)).doit()
diff_g_eq_h

$\displaystyle \frac{\partial}{\partial h} g{\left(h,s \right)} = 20$

diff_g_eq_s = Eq(Derivative(g, s), Derivative(20*h + 2000*s, s)).doit()
diff_g_eq_s

$\displaystyle \frac{\partial}{\partial s} g{\left(h,s \right)} = 2000$

diff_lm_h = Eq(Derivative(100*h**Rational(2, 3)*s**Rational(1, 3), h), lamda*Derivative(20*h + 2000*s, h)).doit()
diff_lm_h

$\displaystyle \frac{200 \sqrt[3]{s}}{3 \sqrt[3]{h}} = 20 \lambda$

diff_lm_s = Eq(Derivative(100*h**Rational(2, 3)*s**Rational(1, 3), s), lamda*Derivative(20*h + 2000*s, s)).doit()
diff_lm_s

$\displaystyle \frac{100 h^{\frac{2}{3}}}{3 s^{\frac{2}{3}}} = 2000 \lambda$

budget_eq = Eq(20*h + 2000*s, B)
budget_eq

$\displaystyle 20 h + 2000 s = 20000$

solns = solve([diff_lm_h, diff_lm_s, budget_eq], [h, s, lamda], dict = True)[0]

R_eq.subs(solns).n()

$\displaystyle R{\left(\frac{2000}{3},\frac{10}{3} \right)} = 11399.8396445113$