Most suitable metric for the Solar system? 
*

*If I wanted to solve the Einstein equations for the solar system, which choice of $g_{\mu\nu}$  and $T_{\mu\nu}$ is more suitable? 

*I thought about using a Schwarzschild metric near each planet, but how to connect them?
 A: Solving the Einstein equation for a system as complex as the Solar System could only be done numerically, and in any case it's not terribly useful. Nothing in the Solar System, is relativistic enough to need more than a linearised treatment (this is how Einstein calculated the precession of Mercury).
Actually even solving Newton's equation for a system as complex as the Solar System can only be done numerically. The way you actually do things like detecting the presence of Neptune is done using perturbation theory. If you were interested in relativistic effects you start with a symmetric solution, apply classical perturbations using Newton's law then finally apply corrections using linearised GR. 
To take the example of Mercury that I mentioned above: if the Solar System consisted only of the Sun and Mercury (and both were perfect spheres) the orbit of Mercury wouldn't precess. However it's observed to precess by 574 arc-seconds per century. 531 arc-seconds of this are down to classical perturbations by other planets, and this was known before GR was formulated. Only 43 arc-seconds was left for Einstein to explain, which he did using a linearised approximation.
A: You just add the metric in isotropic form. You need the isotropic Schwartschild solution, which is found by transforming the r coordinate so that the metric looks like:
$$ - g_{00}(u) dt^2 + g_{uu}(u) (du^2 + u^2 d\Omega^2)$$
This requires choosing $u(r)$ appropriately. Then you note that the result can be transformed again:
$$ - g_{00}(u) dt^2 + g_{uu} (dx^2 + dy^2 + dz^2) $$
Where $u = \sqrt{x^2 + y^2 + z^2}$. Then you write the linearized approximation, which turns out to be:
$$ - dt^2 + dx^2 + dy^2 + dz^2 - h(u) dt^2 - g(u)( dx^2 + dy^2 + dz^2) $$
And then you superpose the solution for the invidual masses, since the linearized small perturbation metric is additive among separate sources (like any other linear field theory). The result can be simply expressed in terms of the Newtonian potential:
$$ h_{tt} = 2 \phi(x,y,z)$$
$$ h_{ii} = 2 \phi(x,y,z) $$
Where $\phi(u)$ is the ordinary Newtonian potential. The superposition is valid in the weak field limit, which for the Solar system is essentially exact. The result can also be derived without going through isotropic coordinates, because the difference in r and u is higher order in the gravitational field, but I prefer to do it this way.
For moving objects, you have to tensorially boost the Schwarzschild metric to a moving frame. This introduces the Gravitomagnetic forces, and it is important in the solar system. The boosting is by the tensor law using the Lorentz boost. This gives the off diagonal h's for a moving source, and it is fine as long as the retardation effects are negligible, which is true throughout the solar system.
This approximation is introduced by Einstein in 1916 in the original papers on GR, and he uses it to calculate the effects of GR in the solar system, a few months before the exact Schwarzschild solution appears. Einstein thought the field equations were too complicated for an exact solution until Schwarzschild proved him wrong.
