You may be interested in the Bertotti-Robinson solution to the field equations. It is the only conformally flat solution to the electrovacuum equations with a non-null (in the sense of the electromagnetic invariants) Faraday tensor. (See also Chapter 7 of Griffiths and Podolsky, Exact space-times in Einstein's general relativity.) (Note that in the form presented in the ScienceWorld link it is not obviously static. See below.)
This solution, however, has a uniform electromagnetic field, and is not asymptotically flat. In the asymptotically flat case, using the Jebsen-Birkhoff-Eiseland rigidity theorem, a spherically symmetric solution (as you postulated in the first line of your question) must admit an additional symmetry, and the field equations reduce to ODEs. Asymptotic flatness allows you to prescribe asymptotic data "at infinity" and integrate back. In this case you have a stronger result that for any (possibly nonlinear) model of electromagnetism that reduces to Maxwell theory in the weak field limit, a spherically symmetric solution must be static and must also be singular along the symmetry axis. For a review of this fact and other related results, see the recent paper by Shadi Tavildar-Zadeh.
Similar to the Bertotti-Robinson solution is also the Melvin Universe, which is a static non-singular, cylindrical symmetric space-time with an electric field directed along the axis of the cylinder. (You can find more about this solution in Griffiths and Podolsky's book mentioned above.) This solution isn't really physical: it was shown by Thorne, and by Melvin and Wallingford that the geometry is very strange in the Melvin Universe: time-like geodesics with only radial components to the velocity cannot escape to radial infinity ever; the geodesics oscillate and repeatedly cross the central symmetry axis. In otherwords, the gravitational effects induced by the electric field in this solution is very strong. (For a non-linear electrodynamics version of the Melvin Universe, there's a paper of Gibbons and Herdeiro on the subject.)
Since you also asked how one coms up with these solutions: a good method is guess and check. Guess an ansatz (for example, a symmetry class), reduce the equation, and try to solve. For example, this would be how I derive the Bertotti-Robinson solution:
We will make the horribly unphysical ansatz that the solution is spherically symmetric and that the space-time is a direct product space-time. So in particular the geometry decomposes into $\mathbb{S}^2 \times (\mathbb{R}^2,g)$. Jebsen-Birkhoff-Eiseland still applies, so we can take one coordinate of $\mathbb{R}^2$ to be the Killing coordinate.
The metric then takes the canonical form
$$ -A(r) dt^2 + B(r) dr^2 + ds^2_{\mathbb{S}^2}$$
As a direct product metric, the Ricci curvatures are just the sums of the component Ricci curvatures. So using the standard formulae for curvatures of two dimensional surfaces, we quickly compute that in the spherical directions the Ricci curvature is that of the standard sphere, which is equal to the metric (2 dimensional manifold with constant curvature 1). In the $\mathbb{R}^2$ direction, the Gaussian curvature is
$$ - \frac{1}{2\sqrt{A(r)B(r)}} \left( \frac{A'(r)}{\sqrt{A(r)B(r)}}\right)' $$
The trace-free condition from Einstein-Maxwell requires that the Gaussian curvature equal to -1, which you can either solve explicitly for by integrating in $r$, or by "knowledge" that the only 2 dimensional Lorentzian manifold with constant negative curvature is AdS2.
So if you take $r$ to be the "inverse radial" coordinate (where $r = 0$ corresponds to spatial infinity), the usual presentation of AdS2 gives that $A(r) = B(r) = r^{-2}$.
For the Maxwell field, you observe that by spherical symmetry $F$ only has $dr\wedge dt$ component (if you toss out magnetic monopoles; you can gain back the other term using a duality rotation anyway). So $F$ is given by a scalar function $f$. The curvature condition implies that $R_{rr} = T_{rr} = F_{rt}F_{rt}g^{tt} = -f^2 A^{-1} = -1g_{rr} = -1/r^2$ (observe that you can say the same thing using $R_{tt}$ or the spherical components). This means that $f^2 = r^{-4}$, or $f = r^{-2}$.
In remains to check that such a solution verifies Maxwell's equation. $dF = f' dr\wedge dr\wedge dt = 0$ clearly. while the Hodge dual $*F = \sqrt{|g|}g^{tt}g^{rr} f dvol_{\mathbb{S}^2} = - dvol_{\mathbb{S}^2}$, so $d*F = 0$ also. So the Bertotti-Robinson solution is indeed a static solution to the field equations. That the solution is non-singular follows from the fact that it is a direct product of two constant curvature spaces, so the curvature is bounded, and that the electromagnetic invariant $g(F,F)$ is constant valued, so the Maxwell field does not blowup either.
:-)(And if you allow cosmological constants, so are dS and AdS...) – Willie Wong Jun 15 '11 at 16:36