There is only one Feynman diagram responsible for $e^+ e^- \to \mu^+ \mu^-$ process in the leading order. The following Redberry code produces squared matrix element of this process:
//setting up matrices
//gamma, vertex
defineMatrices 'G_a', 'V_i', Matrix1.matrix,
//electron & muon wave functuins
'v[p_a]', 'u[p_a]', Matrix1.vector,
//their conjugations
'cv[p_a]', 'cu[p_a]', Matrix1.covector
//photon propagator
def G = 'G_mn[k_a] = -I*g_mn/(k_a*k^a)'.t
//vertex
def V = 'V_i = -I*e*G_i'.t
//matrix element
def M = 'cv[p2_a]*V_i*u[p1_a]*G^ij[p1_a + p2_a]*cu[k1_a]*V_j*v[k2_a]'.t
//substitute Feynman rules
M = (V & G) >> M
//list of Mandelstam & mass shell substitutions
def mandelstam = setMandelstam(
['p1_m': 'me', 'p2_m': 'me', 'k1_m': 'mu', 'k2_m': 'mu'])
//simplify matrix element
M = (EliminateMetrics & ExpandDenominator & mandelstam) >> M
//complex conjugation of matrix element
def MC = Conjugate >> M
MC = 'u[p1_a]*cv[p2_a] = v[p2_a]*cu[p1_a]'.t >> MC
MC = 'v[k2_a]*cu[k1_a] = u[k1_a]*cv[k2_a]'.t >> MC
//squared matrix element
def M2 = ExpandAll >> (M * MC / 4)
//sum over electron and muon polarizations
M2 = 'u[p1_a]*cu[p1_a] = me + p1_a*G^a'.t >> M2
M2 = 'u[k1_a]*cu[k1_a] = mu + k1_a*G^a'.t >> M2
M2 = 'v[p2_a]*cv[p2_a] = -me + p2_a*G^a'.t >> M2
M2 = 'v[k2_a]*cv[k2_a] = -mu + k2_a*G^a'.t >> M2
//trace of gamma matrices
M2 = DiracTrace['G_a'] >> M2
//simplifications
M2 = (ExpandAndEliminate & mandelstam) >> M2
M2 = 'u = 2*(mu**2 + me**2) - s - t'.t >> M2
M2 = Factor >> M2
println M2
> 2*(2*mu**4+2*t**2+2*s*t+2*me**4-4*t*mu**2-4*me**2*t+4*me**2*mu**2+s**2)*
s**(-2)*e**4
Which is a well known formula:
\[
\frac{2 e^{4}}{s^2}\, (4 m_u^{2} m_e^{2}+2 t^{2}+2 m_u^{4}+s^{2}+2 m_e^{4}-4 m_e^{2} t-4 m_u^{2} t+2 s t)
\]