Does this polynomial evaluate to prime number whenever $x$ is a natural number?

I am trying to prove or disprove following statment:

$x^2-31x+257$ evaluates to a prime number whenever $x$ is a natural number.

First of all, I realized that we can't factorize this polynomial using its square root like $$ax^2+bx+c=a(x-x_1)(x-x_2)$$ because discriminant is negative, also I used excel for numbers from 1 to 1530 (just for checking), and yes it gives me prime numbers, unfortunately I dont know what the formula is, which evaluates every natural number for given input, maybe we may write $x_{k+1}=k+1$ for $k=0,\ldots\infty$, but how can I use this recursive statment? I have tried instead of $x$, put $k+1$ in this polynomial, and so I got

$$(k+1)^2-31k+257=k^2+2k+1-31k+257=k^2-29k+258$$ but now this last polynomial for $k=1$ evaluates to $259-29=230$, which is not prime, but are my steps correct? Please help me.

If $x$ is divisible by $257$ then so is $x^2 - 31 x + 257$. More generally, if $f(x)$ is any polynomial $f(x)$ is divisible by $f(y)$ whenever $x-y$ is divisible by $f(y)$. So there are no non-constant polynomials that produce primes for all positive integers.

There are obvious counterexamples [hint: $\,f(0)=p\Rightarrow p\mid f(pn)\,$] and this idea leads to a general proof that no such prime-producing polynomial exists.

However, you probably wonder why the first $32$ values are all prime. To explain that, notice that if you shift x by $16$ then it becomes $\rm\:f(x) = x^2 + x + 17.\:$ Now one may apply the following theorem on prime-producing polynomials (generalizing the famous example $\rm\:x^2+x+41\:$ noticed by Euler), noting that the ring of integers of $\rm\:\mathbb Q(\sqrt{-67})\:$ is a UFD.

Theorem $\ $ The polynomial $\rm\ f(x)\ =\ (x-\alpha)\:(x-\alpha')\ =\ x^2 + x + k\ $ assumes only prime values for $\rm\ 0\ \le\ x\ \le\ k-2 \ \iff\ \mathbb Z[\alpha]\ $ is a PID.

Hint $\: (\Rightarrow)\ $ Show all primes $\rm\ p\le\sqrt{n},\, n=1\!-\!4k\,$ have $\rm \left[\frac{n}p\right]\! = -1,\,$ i.e. $\rm\! \bmod p\!:\rm\, n\,$ is nonsquare, so no primes split/ramify.

For proofs, see e.g. Cohn, Advanced Number Theory, pp. 155-156, or Ribenboim, My numbers, my friends, 5.7 p.108. Note: both proofs employ the bound $\rm\ p < \sqrt{n}\ $ without explicitly mentiioning that this is a consequence of the Minkowski bound - presumably assuming that is obvious to the reader based upon earlier results. Thus you'll need to read the prior sections on the Minkowski bound. Compare Stewart and Tall, Algebraic number theory and FLT, 3ed, Theorem 10.4 p.176 where the use of the Minkowski bound is mentioned explicitly.

Alternatively see the self-contained paper [1] which proceeds a bit more simply, employing Dirichlet approximation to obtain a generalization of the Euclidean algorithm (the Dedekind-Rabinowitsch-Hasse criterion for a PID). If memory serves correct, this is close to the approach originally employed by Rabinowitsch when he first published this theorem in 1913.

[1] Daniel Fendel, Prime-Producing Polynomials and Principal Ideal Domains,
Mathematics Magazine, Vol. 58, 4, 1985, 204-210

For $x=32$, we have $$x^2-31x+257=289=17^2.$$

An easy way with polynomials in a single variable is to set the variable $x$ equal to the constant term - in this case $x=257$. If this turns out negative*, choose a suitably high multiple of the constant term (this assumes the leading coefficient is positive). The factorisation is obvious and you don't need to multiply large numbers.

*or equal to the constant when the constant is prime