The classical proof via differential geometry goes like this:  

Suppose that the surface in $\mathbb{R}^3$ is smooth and parametrize it locally in the form $X(s,t)$ where the two rulings are defined by holding either $s$ or $t$ constant.  This is a local argument, so, for simplicity, I'll assume that the domain of $X$ is a rectangle in the $st$-plane.  Of course, one can reparametrize in $s$ and/or $t$ separately, and this will turn out to be useful at some point in the calculation. 

The two tangent vector fields $X_s$ and $X_t$ are linearly independent and are the tangents to the two rulings.  Since $X_{ss}$ is the acceleration of the $t$-ruling, it follows that $X_{ss} = f X_s$ for some function $f$.  Similarly $X_{tt} = g X_t$ for some function $g$.
Note that, if one reparametrized, using $\bar s$ and $\bar t$ instead of $s$ and $t$, then one would have
$$
X_{\bar s} = \frac{ds}{d\bar s}\ X_s\quad\text{and}\quad
X_{\bar s\bar s} = \left(f + \frac{d^2s}{d\bar s^2}\left(\frac{ds}{d\bar s}\right)^{-1}\right)\ X_{\bar s},
$$
with similar formulae for $X_{\bar t}$ and $X_{\bar t\bar t}$.  This will be useful below.


Since the surface does not lie in a plane, $X_{st}$ cannot be a linear combination of $X_s$ and $X_t$ (otherwise, the plane spanned by $X_s$ and $X_t$ would be fixed, and the surface would lie in plane). This means that $X_s$, $X_t$, $X_{st}$ is a basis of $\mathbb{R}^3$, and, as such, there are equations of the form
$$
\begin{pmatrix} dX_s& dX_t & dX_{st}
\end{pmatrix}
= \begin{pmatrix} X_s& X_t & X_{st}
\end{pmatrix}
\begin{pmatrix} f\ ds& 0 & f_t\ ds\\\\
0 & g\ dt & g_s\ dt\\\\
dt & ds & f\ ds + g\ dt
\end{pmatrix}
$$
(The equations for $dX_{st}$ follow since $(X_{st})_s = (f X_s)_t = f_t\ X_s + f\ X_{st}$, etc.)  By comparing partials, or by using the structure equations above (i.e., expanding out the consequences of $d(d(X_{st}))=0$, etc.), one sees that $d(f\ ds + g\ dt) = 0$, so that there must exist a function $h$ such that $f = 2h_s$ and $g = 2h_t$.  (The coefficient of $2$ avoids some fractions later.)
The equation now becomes
$$
\begin{pmatrix} dX_s& dX_t & dX_{st}
\end{pmatrix}
= \begin{pmatrix} X_s& X_t & X_{st}
\end{pmatrix}
\begin{pmatrix} 2h_s\ ds& 0 & 2h_{st}\ ds\\\\
0 & 2h_t\ dt & 2h_{st}\ dt\\\\
dt & ds & 2h_s\ ds + 2h_t\ dt
\end{pmatrix}
$$
Moreover, the structure equations now imply that $d(e^{-2h}h_{st})=0$, so $h_{st} = C\ e^{2h}$ for some constant $C$.  By adding a constant to $h$, one can reduce to the case that $C$ is one of $0$, $1$, or $-1$.

Consider the case $C=0$ (which needs to be treated separately in any case).  Then $h_{st}=0$, so that, in particular, $f=2h_s$ is a function of $s$ alone and $g=2h_t$ is a function of $t$ alone. Using the change of variables formulae mentioned above, one can then change variables in $s$ so as to arrange that $f = 0$ and change variables in $t$ to arrange that $g=0$.  Thus, the equations have reduced to
$$
d(X_s) = X_{st}\ dt,\qquad d(X_t) = X_{st}\ ds,\qquad d(X_{st})=0.
$$
Thus $X_{st} = v_3$ where $v_3$ is a constant vector.  Then $d(X_s-tv_3) = 0$ and $d(X_t - s v_3) = 0$, so there exist constant vectors $v_1$ and $v_2$ so that $X_s = v_1 + t v_3$ and so that $X_t = v_2 + s v_3$.  Finally, this implies that 
$$
dX = X_s\ ds + X_t\ dt = (v_1+tv_3)\ ds + (v_2 + s v_3)\ dt 
= d\bigl(sv_1+t v_2 + st v_3 \bigr), 
$$
so that
$$
X = v_0 + sv_1+t v_2 + st v_3
$$
for some constant vector $v_0$.  Thus, $X(s,t)$ parametrizes a hyperbolic paraboloid.

Now, consider the case $C\not=0$.  The structure equations have become

$$
\begin{pmatrix} dX_s& dX_t & dX_{st}
\end{pmatrix}
= \begin{pmatrix} X_s& X_t & X_{st}
\end{pmatrix}
\begin{pmatrix} 2h_s\ ds& 0 & 2C\ e^{2h}\ ds\\\\
0 & 2h_t\ dt & 2C\ e^{2h}\ dt\\\\
dt & ds & 2\ dh
\end{pmatrix}
$$
where $h_{st} =  C\ e^{2h}$.  Moreover, one calculates
$$
d(e^{-2h} X_{st}) = 2C\ (X_s\ ds + X_t\ dt) = 2C\ dX
$$,
so it follows that $e^{-2h}X_{st} = 2C X + v_0$ for some constant vector $v_0$.  In particular, showing that the vector-valued function $E_3 = e^{-2h}X_{st}$
takes values in a hyperboloid of $1$-sheet will finish the proof.  

To this end, consider the new frame field
$$
\begin{pmatrix}E_1 & E_2 & E_3\end{pmatrix}
= \begin{pmatrix}e^{-h}X_{s} & e^{-h}X_{t} & e^{-2h}X_{st}\end{pmatrix}.
$$
Calculation shows that it satisfies the structure equation
$$
\begin{pmatrix} dE_1& dE_2 & dE_3
\end{pmatrix}
= \begin{pmatrix} E_1& E_2 & E_3
\end{pmatrix}
\begin{pmatrix} h_s\ ds-h_t\ dt& 0 & 2C\ e^h\ ds\\\\
0 & h_t\ dt - h_s\ ds & 2C\ e^h\ dt\\\\
e^h\ dt & e^h\ ds & 0
\end{pmatrix}.
$$
Note that the $3$-by-$3$ matrix of $1$-forms on the right takes values in the vector space
$\frak{g}$ consisting of matrices of the form
$$
\begin{pmatrix} x_1& 0 & 2C\ x_2\\\\
0 & -x_1 & 2C\ x_3\\\\
x_3 & x_2 & 0
\end{pmatrix}.
$$ 
This is, of course, the Lie algebra of the subgroup $O(Q)\subset GL(3,\mathbb{R})$ consisting of the matrices that satisfy $A^TQA = Q$ where
$$
Q = \begin{pmatrix} 0& 1 & 0\\\\
1 & 0 & 0\\\\
0 & 0 & -2C
\end{pmatrix}
$$
It follows that there is an invertible linear transformation $L$ of $\mathbb{R}^3$ such that the matrix $LE$ takes values in $O(Q)$, where $E = (E_1\ E_2\ E_3)$.  In particular $L$ carries the image of $E_3$ into the hyperboloid of $1$-sheet $2x_1x_2 - 2C x_3^2 = -2C$.  By the remark above, it follows that $X(s,t)$ must be the image of this quadric under an affine transformation of $\mathbb{R}^3$, as was to be shown.