Suppose that I am given a subscheme $Y$ of $\mathbf{P}^n_{\mathbf{Z}}$, flat over $\operatorname{Spec}\mathbf{Z}$ and with smooth generic fiber $Y_{\mathbf{Q}}$, defined by the vanishing of some homogeneous polynomials $$ F_1, \ldots, F_k \in \mathbf{Z}[X_0,\ldots,X_n]. $$ How does one determine the set $S$ consisting of primes $p$ such that the fiber $Y_{\mathbf{F}_p}$ is non-smooth?

The way I'm going about the problem seems a bit crude to me, and has some problems. In particular, I can't prove that my algorithm solves the problem.

For simplicity suppose that $k=1$, and write $F=F_1$. In a nutshell, I compute the derivatives $f_i := \frac{\partial F}{\partial X_i}$, and determine an integer $N$ such that $N$ is contained in the ideal $I = (F, f_1, \ldots, f_n) \subset \mathbf{Z}[X_0,\ldots,X_n]$ by repeatedly taking resultants of polynomials in $I$. Then at least for all $p$ not dividing $N$, the fiber $Y_{\mathbf{F}_p}$ is smooth, since for such $p$ the ideal $I_p \subset \mathbf{F}_p[X_0,\ldots,X_n]$ generated by the reduced polynomials $\widetilde{F}$, $\widetilde{f}_1$,$\ldots$,$\widetilde{f}_n$ contains $1$.

One problem is that it's quite expensive computationally. I'm using resultants to eliminate the variables $X_i$ one by one, say in the order $i=0,1,\ldots,n$. Then every time I've eliminated one of the $X_i$s, I get a lot of polynomials that don't contain $X_0$ up to $X_i$, and I take them all into the next elimination round. (You have to be careful when the resultant comes out $0$, and perhaps there are some other subtleties that I'm forgetting, but this is basically it.) Another problem with this is that the final outcome seems to depend on the order in which I take the $i$s. (This may have something to do with singularities being located at infinity with respect to one of the $X_i$s, but I really can't see the geometry of what I'm doing clearly enough to be confident about this.) And I don't see why my method necessarily gives me the whole $S$ - instead of some upper bound on it - even when all orderings of the $X_i$ are taken into account.

My questions:

1. Does my method give the whole set $S$?
2. If yes, why? If no, what method does work?
3. Regardless of the answers to 1. and 2., is there a faster or more natural way of finding the primes of bad reduction?

Suppose that I am given a subscheme $Y$ of $\mathbf{P}^n_{\mathbf{Z}}$, flat over $\operatorname{Spec}\mathbf{Z}$ and with smooth generic fiber $Y_{\mathbf{Q}}$, defined by the vanishing of some homogeneous polynomials $$ F_1, \ldots, F_k \in \mathbf{Z}[X_0,\ldots,X_n]. $$ How does one determine the set $S$ consisting of primes $p$ such that the fiber $Y_{\mathbf{F}_p}$ is non-smooth?

The way I'm going about the problem seems a bit crude to me, and has some problems. In particular, I can't prove that my algorithm solves the problem.

For simplicity suppose that $k=1$, and write $F=F_1$. In a nutshell, I compute the derivatives $f_i := \frac{\partial F}{\partial X_i}$, and determine an integer $N$ such that $N$ is contained in the ideal $I = (F, f_1, \ldots, f_n) \subset \mathbf{Z}[X_0,\ldots,X_n]$ by repeatedly taking resultants of polynomials in $I$. Then at least for all $p$ not dividing $N$, the fiber $Y_{\mathbf{F}_p}$ is smooth, since for such $p$ the ideal $I_p \subset \mathbf{F}_p[X_0,\ldots,X_n]$ generated by the reduced polynomials $\widetilde{F}$, $\widetilde{f}_1$,$\ldots$,$\widetilde{f}_n$ contains $1$.

One problem is that it's quite expensive computationally. I'm using resultants to eliminate the variables $X_i$ one by one, say in the order $i=0,1,\ldots,n$. Then every time I've eliminated one of the $X_i$s, I get a lot of polynomials that don't contain $X_0$ up to $X_i$, and I take them all into the next elimination round. (You have to be careful when the resultant comes out $0$, and perhaps there are some other subtleties that I'm forgetting, but this is basically it.) Another problem with this is that the final outcome seems to depend on the order in which I take the $i$s. (This may have something to do with singularities being located at infinity with respect to one of the $X_i$s [edit: on second thought, this doesn't make any sense; see comment by François], but I really can't see the geometry of what I'm doing clearly enough to be confident about this.) And I don't see why my method necessarily gives me the whole $S$ - instead of some upper bound on it - even when all orderings of the $X_i$ are taken into account.

My questions:

1. Does my method give the whole set $S$?
2. If yes, why? If no, what method does work?
3. Regardless of the answers to 1. and 2., is there a faster or more natural way of finding the primes of bad reduction?