Why does this manual bilinear transformation yield different results from Matlab's?
A couple things:

Before making the substitution $s = \frac{2}{T} \frac{z-1}{z+1}$, you need to prewarp the cutoff frequency by making the substitution:

$$
\omega_{c,w} = \frac{2}{T} \tan (\omega_c\frac{T}{2})
$$

where $\omega_{c,w}$ is the warped cutoff frequency. This is necessary because the bilinear transform maps the left-half plane in the Laplace domain (used in analog filter design) to the unit circle in the $z$-domain in a nonlinear fashion. Therefore, as you approach the Nyquist rate (digital frequencies of $\pm \pi$), the approximation to the analog filter prototype becomes inaccurate.

Also, the second parameter that you are passing to the butter function is the normalized cutoff frequency, not the sampling interval $T$. The normalized frequency used by that function is in the interval $(0,1)$ and is equal to the ratio of the desired cutoff frequency to the Nyquist rate:

$$
\omega_n = \frac{\omega_c}{2 \pi \frac{f_s}{2}}
$$

$$
\omega_n = \frac{\omega_c}{\pi f_s}
$$

$$
\omega_n = \frac{\omega_c T}{\pi}
$$