Evaluating a series of Legendre polynomials
Compute legendre(n) once per loop iteration instead of twice.
Use np.hypot instead of np.sqrt(x**2 + z**2).
Don't use phivec, just pass the arrays X and Y directly to phi, which is already almost vectorized. The only change that's needed is to use np.minimum and np.maximum instead of min and max respectively.


Revised code:

def phi(x, z, nlimit=35):
    r = np.hypot(x, z)
    costh = z / r
    rmin = np.minimum(r, 1)
    rmax = np.maximum(r, 1)
    dr = rmin / rmax
    suma = 0.0
    for n in range(nlimit):
        Ln = sp.legendre(n)
        suma += dr**n * Ln(costh) * Ln(0.0)
    return suma / rmax

X, Y = np.meshgrid(*[np.linspace(-2, 2, 50)] * 2)
Z = phi(X, Y)


I find that this takes about 0.05 seconds, which ought to be fast enough for you.