How do i solve a 2 variable polynomial over 1 variable

edit

The polynomial Phi can only be evaluated at two inputs (using the function call syntax). The method roots is only available for single variable polynomials. Your Phi(33,Y) is (in SageMath) still a multivariable polynomial (in which X does not appear). To convert it to a single variable polynomial in Y, use the univariate_polynomial method:

sage: f = Phi.subs(X=33).univariate_polynomial(); f
Y^3 + 6150*Y^2 + 5541*Y + 1175
sage: f.parent()
Univariate Polynomial Ring in Y over Finite Field of size 8009
sage: f.roots()
[(898, 1)]
sage: f.roots(multiplicities=False)
[898]

Alternatively:

sage: R.ideal([Phi, X-33]).variety()
[{Y: 898, X: 33}]

This is because Phi(33, Y) is still a polynomial in two variables:

sage: Phi(33, Y).parent()
Multivariate Polynomial Ring in X, Y over Finite Field of size 8009

So, you have to turn it into a one-variable polynomial first:

sage: S.<Y> = GF(8009)[]
sage: S(Phi(33, Y))
Y^3 + 6150*Y^2 + 5541*Y + 1175

sage: S(Phi(33, Y)).parent()
Univariate Polynomial Ring in Y over Finite Field of size 8009

sage: S(Phi(33, Y)).roots()
[(898, 1)]