I'm going to give the same answer as Michael, but phrased somewhat differently: I'll work in a coordinate system whose $xy$-plane is the plane of the orbit of the earth, and whose $z$-axis is the perpendicular to this plane. The summer solstice occurs when the earth's center has largest possible $x$-coordinate, and the polar axis of the earth is therefore tilted so that at the summer solstice, it's in the xz-plane, with the north pole nearer to the sun, and the south pole farther from the sun.
Letting $\theta = \dfrac{2\pi(\text{date} - \text{summerSolstice})}{365}$, we have the position of the center of the earth as
$$
C(\theta) = (R \cos \theta, 0, R \sin \theta)
$$
where $R$ is the radius of the earth's orbit, which will not be needed further.
$$
\newcommand{\vv}{\mathbf v}
\newcommand{\uu}{\mathbf u}
\newcommand{\ww}{\mathbf u}
\newcommand{\ey}{\mathbf {e_2}}
$$
The vector along the earth's axis, from S to N, is
$$
\vv = (-\sin \alpha, 0, \cos \alpha),
$$
where $\alpha = 23.56 \text{ degrees}$ is the inclination of the axis. (You can make $\alpha$ a function of time to include nutation, etc., but I'm not going to do that). The vector perpendicular to this, in the $xz$ plane, pointing mostly away from the sun at the summer solstice, is
$$
\uu = (\cos \alpha, 0, \sin \alpha).
$$
A third vector, perpendicular to both of these, is the $y$-axis, $\ey$. A typical point $P$ on the earth can be described by its displacement from the earth's center, i.e., by
$$
\ww = P - C(\theta)
$$
If the latitude of $P$, measured between $\pi/2$ at the south pole to $\pi/2$ at the north pole, is $\lambda$, then the vector $\ww$ looks like
$$
\ww = r\sin (\lambda) [ \sin(t) \ey + \cos(t) \uu] + r\cos (\lambda) \vv
$$
where $t$ is the time of day, measured from $0$ to $2 \pi$, and $r \ll R$ is the radius of the earth.
Sunrise/sunset occur when the vector $\ww$ is orthogonal to the ray from the origin to the point $P$, i.e., when $\ww$ is orthgonal to $\ww + C(\theta)$, i.e.
$$
\ww \cdot \ww + \ww \cdot C(\theta) = 0.
$$
At this point, there are two quite reasonable approximations to be made.
Pretend that the time $t$ and the date $\theta$ are independent. We fix $\theta$ at the angle for mid-day, on the grounds that it's changing very slowly.
The dot product of $\ww$ with itself is small compared to the dot product of $C(\theta)$ with $\ww$, since the magnitude of $C(\theta)$ is the radius of the earth's orbit, while the magnitude of $\ww$ is the radius of the earth, so pretend that the $\uu \cdot \uu$ term is missing in the equation.
The first of these (and perhaps the second) is used in Michael Hardy's solution. As I recall from my study of celestial navigation, the impact is fairly small, although it does vary from day to day once we no longer assume a circular orbit, so it's not to be ignored if you want ultimate accuracy. The second introduces an error in angle no larger than $r/R$, which is something like $1/2000$ radian, although this is a coarse overestimate; that's about 1.7 minutes of arc, which would turn into about 6 minutes of time.
At any rate, one can write $t$ in terms of $\theta$: multiply $\theta$ by 365; divide by $2 \pi$, and then take the result mod 1, and multiply by $2 \pi$ to get $t$.
Once you've done that, the equation above involves only $\theta$, and you can solve for when it's zero (although you probably have to resort to numerical means; exact trig isn't likely to work out nicely).
The only real difference of this from the previous solution is that time and date get related, and the computation is done with dot products rather than classical geometry.