Such a simple question, but it opens up so many cans of worms. Here's my crack at a comprehensive answer simpler than what you'd find in a textbook. I'm sure a numerical relativist would have a much more inside scoop... this is coming from more of a mathematical perspective. First I'll try to clear up some general difficulties which I think were misunderstood in the question, then get to the example at the end. I was excited to see this question especially since, as you said, there's a big lack of simple examples here. Had a little too much fun, I hope the answer's not too long to be helpful!
Basically, yes, what you want to do can be done. But if you want to include a nontrivial matter distribution, it gets really hard really fast. To see why, we have to look at the overall problem we're trying to solve. It's not just "solve Einstein's field equations". It's "solve Einstein's field equations coupled to matter equations". We're trying to minimize the action
$$ S = S_{gravity} + S_{matter}. $$
Varying the action with respect to the metric gives Einstein's equation, while varying with respect to whatever variables are in $S_{matter}$ gives the matter time evolution equations (this might be the curved-spacetime Klein-Gordon equation for a scalar field, for example).
Let's suppose, for definiteness, that the matter content is going to be modeled by some classical fields. Varying the total action gives a second-order coupled system of nonlinear PDE's, some coming from Einstein's equation, some coming from matter equations. Since the system is second order, we basically need to provide an initial configuration and initial time derivative for both the metric and each of the matter fields.
So even before trying to actually solve the problem, it's not as simple as: "Let me just specify $T_{\mu\nu}$ and then solve Einstein's equation for the metric". First you have to choose how to model the matter, and hope it's going to evolve the way you want it to. Then you have to set the initial data for the metric, and choose initial data for the matter fields that give you the desired initial $T_{\mu\nu}$. There are probably multiple ambiguous ways to do this, so watch out. Now you have to check that all your initial data match any constraints on the system (e.g. gauge fixing constraints). If you managed all that, you can finally hit go and pray that the time evolution does something like what you expected. Clearly this is a pretty big undertaking, and we haven't even considered all the hairy details yet. So much for simple examples.
And, by the way, you can't just ignore the matter dynamics and treat $T_{\mu\nu}$ as a fixed source, since, as Elio Fabri has pointed out in his answer, you can't really specify $T_{\mu\nu}$ without knowing the spacetime structure.
Fortunately, however, you asked about the Schwarzschild solution, and the Schwarzschild solution doesn't have any matter in it! Well, actually, the matter is all relegated to "elsewhere", and manifests as a boundary condition on the Weyl tensor. But anyway, locally, no matter. That sounds a lot easier. So let's talk about vacuum solutions, where $T_{\mu\nu}=0$ and $S=S_{gravity}$ only.
Vacuum Solutions
First of all, we need to pose the question.
Let's work in $D=N+1$ spacetime dimensions.
We are given:
- An initial N-dimensional manifold $\Sigma$.
- An initial spatial (positive-definite) metric $h_{ij}$ on $\Sigma$.
- An initial time derivative of the metric, usually specified by the extrinsic curvature tensor $K_{ij}$ on $\Sigma$.
- Some kind of spatial boundary conditions.
We want to obtain a spacetime $M$ with metric $g_{\mu\nu}$ such that:
- The spacetime obeys the vacuum Einstein equation $G_{\mu\nu}=0$.
- You can slice up $M$ into spacelike hypersurfaces, each topologically equivalent to $\Sigma$, and each obeying the desired boundary conditions, such that one of the slices is equivalent to the initial data.
- The spatial hypersurfaces foliating $M$ can be labeled smoothly by a parameter $t$, increasing toward the future.
To convert this into a differential equation for the spatial metric evolving on $\Sigma$ we have to write out
- $G_{\mu\nu}=0$,
- The "Gauss-Codazzi" constraints relating $h_{ij}$ to $K_{ij}$ on $\Sigma$ (these come from differential geometry of embedded surfaces),
- Four extra (gauge) constraints accounting for the coordinate ambiguity,
in terms of $h_{ij}$, $K_{ij}$, and a parameter $\alpha$ called the lapse function (which is basically just $g_{00}$, and tells you how much proper time elapses while traveling normal to $\Sigma_t$).
After all that, you finally have a set of differential equations for the spatial metric, along with some constraints on the initial data.
Yikes, still got complicated pretty fast, huh?
If you want to do all this without a giant effort, it's only going to be in the absolute simplest cases, and you basically have to already know the answer. We'll take a look at that, but first we need to consider some more of the challenges.
EDIT: Thought of a better way to explain previous section, but erasing that much would be bad form, so here's a quick addendum to it.
Ansatz
Assume that the solution spacetime to be obtained has the topology $\mathbb{R} \times \Sigma$ and can be covered by coordinates $(t,x^i)$ with metric
$$ ds^2 = - \alpha^2 \, dt^2 + h_{ij} \, (dx^i + \beta^i dt) (dx^j + \beta^j dt). $$
We call $\alpha$ and $\beta^i$ the lapse function and shift vector, respectively.
If we can get a spacetime metric of this form which
- obeys $G_{\mu\nu}=0$,
- matches our initial data on the slice $t=0$,
- matches any other constraints or boundary conditions,
we'd be done.
The equations of motion (see example at the the end for $\beta^i=0$ case) are derived by calculating extrinsic curvature in the assumed metric, applying $G_{\mu\nu}=0$ in the assumed metric, and utilizing the Gauss-Codazzi relations. The timelike components of $G_{\mu\nu}=0$ are not dynamical, and are called the Hamiltonian and Vector constraints.
Solutions to the equations of motion are not unique unless you specify four additional "gauge" constraints on the spacetime metric. The simplest possible choice would be $\alpha=1$ and $\beta^i=0$. But often more complicated constraints, which may involve the metric and extrinsic curvature in addition to the lapse and shift, make the problem easier.
Slicing: A big can of worms
Fundamentally, the question of how spacetime is sliced into spatial parts is what makes the problem so hard.
To see why, let's consider the simplest possible example: flat, empty, Minkowski spacetime.
Instead of starting from scratch, work backwards. We'll start with the full spacetime, slice it up into time-dependent spatial slices, choose one slice to act as the initial data, integrate from that initial data, and see if we recover the initial slicing, and thereby recover the Minkowski spacetime (I'm not going to actually perform the integration now, but that's the idea we're considering).
Consider the following four slicings:
Why is the difference between these so important? Consider each one:
(a) If you slice like this, the spatial metric is always flat, and the integration is easy. You can only do this because you already knew the answer.
(b) If you slice like this, the spatial metric could be something horrible, and the integration could be difficult, even though the spacetime is trivial. You might do this by accident if you start from scratch.
For many reasons, it's usually impossible to have a single numerical calculation cover the "whole" spacetime (in Minkowski it is possible, of course). This might be because you only know initial data in a certain region of space, because the spacetime itself has no global Cauchy surface, or because the spatial or spacetime manifold has no global coordinate system. When you're only dealing with part of the spacetime, the solution is only valid within a diamond-shaped (in causal diagrams) region called the "domain of dependence" of the surface: basically the subregion of spacetime where no outside influences have time to come in. For the next two we look at slicings of the domain of dependence of an initial surface.
(c) Slicing like this could work, and would make the integration easy, but leaves an ambiguity in the boundary conditions. You could do it, but you have to make sure the numerical solution includes the domain of dependence of the initial surface, and discard parts of the solution outside.
(d) This slicing has obvious boundary conditions and a clear physical interpretation in terms of the domain of dependence, but the spatial metric and differential equations are much worse.
When you start from scratch, you might not even know which of these four cases you're dealing with.
More Worms: Specifying the initial manifold
There's a whole bunch of constraints on specifying initial data for the metric. They come from the coordinate gauge constraints, from the "Gauss-Codazzi" relations I mentioned earlier, and from the timelike components of Einstein's equation.
Let's forget about all that. Much more fundamentally, let's talk about just picking a manifold and throwing a metric on it, with the hope of it corresponding to some particular thing. You can't just do it willy-nilly.
You might think, let me just take the manifold $\mathbb{R}^2$, give it coordinates $(x,y)$, and put the metric $ds^2 = dx^2 + \sin^2 \! x \, \, dy^2$ on it everywhere. Nope! That's the metric of a sphere. It curves back on itself. The manifold you specified doesn't actually make sense, because the geometry and the topology of the manifold are mismatched. There's some fancy theorems about this (sphere theorem, Gauss-Bonnet), but in general the point is that you can't just assume a metric on some coordinate patch makes sense. Locally you're always fine, if you work in a small enough neighborhood. But if you extend the patch far enough, you can run into trouble. The same type of thing can sometimes happen when, for example, you try to define "regularized" versions of Schwarzchild spacetimes that pass through $r=0$.
I don't think this is usually a practical issue in numerical calculations, but I've run into it before on occasion and wanted to mention it because it's pretty neat.
Schwarzschild Example
There's pretty much just one way to do the Schwarzschild example that makes any sense to attempt as a simple case. That is, to cover the entire "exterior" Schwarzschild patch, parameterized by $2M < r < \infty$ and $ -\infty < t < \infty$, sliced on lines of $t=const$.
This is a good choice because it exactly covers the domain of dependence of the initial surface $t=0$, it has obvious spatial boundary conditions, and the metric and curvature are simple.
On the other hand, it's not very physical. For the vacuum exterior to either a black hole or star, we'd have to adjust the initial domain of $r$ in a way that would make the numerical calculation less trivial.
Let's go for it. The spacetime metric we're supposed to end up with is
$$ds^2 = -f(r) \, dt^2 + f(r)^{-1} \, dr^2 + r^2 \, d\Omega^2,$$
where $f(r) = 1-2M/r$.
Foliate on slices $t=const$. The spatial metric on each slice is $ds^2 = f(r)^{-1} dr^2 + r^2 d\Omega^2$. The intrinsic Ricci curvature on each slice is $R_{ij} = \frac{2M}{r^3} \, \textrm{diag}(-f(r)^{-1}, \frac{r^2}{2}, \frac{r^2}{2} \sin^2 \theta)$. The extrinsic curvature on each slice is $K_{ij}=0$ (this isn't trivial).
We're integrating to find the spacetime metric [1]
$$ ds^2 = - \alpha^2 \, dt^2 + h_{ij} \, (dx^i + \beta^i dt) (dx^j + \beta^j dt). $$
Call the spatial coordinates $(x^1,x^2,x^3) = (r, \theta, \phi)$.
As our coordinate gauge conditions choose $\alpha = f(r)^{1/2}$ and $\beta^i=0$ (the $\beta^i$ are "shift vectors", which I left out earlier for simplicity).
For spatial boundary conditions, $h_{ij}$ is periodic in $\theta$ and $\phi$ in the usual spherical way, and $h_{ij}$ is equal to the initial value at $r = 2M$ and $r \to \infty$.
The initial data is $h_{ij} = \textrm{diag}(f(r)^{-1},r^2,r^2 \sin^2 \theta)$ and $K_{ij}=0$.
The initial data already automatically obey all other constraints, since we initally got them by slicing a spacetime.
From [1], the equations of motion with the above setup are
$$ \partial_t h_{ij} = -2 \alpha K_{ij},$$
$$ \partial_t K_{ij} = \alpha [ R_{ij} - 2 K_{il} {K^{l}}_{j} + K K_{ij}] - D_i D_j \alpha, $$
where $K = {K^i}_i$, and $D_{i}$ and $R_{ij}$ are the spatial covariant derivative and intrinsic Ricci curvature on the slice.
At the first timestep, the right hand side of both equations of motion are zero (it's a bit tedious and you need to calculate the Christoffel symbols for the spatial metric on the slice, but the two nonzero terms in the bottom equation eventually cancel out). Since nothing changes, at each subsequent timestep the right hand sides remain zero, and the metric and extrinsic curvature remain constant. We can do this for all $t$.
So the answer we get by integrating is the above spacetime metric, with $h_{ij}(t) = h_{ij}(0)$.
That's the right answer! Exterior Schwarzschild is static! We're done! Pretty trivial answer, but that's how it goes. Clearly, a lot of the steps we took were only simple because we already knew the answer.
Thanks for the opportunity to write this out!
Here's the references I consulted:
[1] Luis Lehner. Numerical Relativity: A review. 2001. https://arxiv.org/abs/gr-qc/0106072.
[2] Wald. General Relativity. Chapter 10.
[3] Hawking and Ellis. Large Scale Structure of Spacetime. Beginning of chapter 7 (it gets ridiculous after the nice intro bit in 7.1).
