Difference between MATLAB's $\tt butter$ and my calculated coefficient

Question

So I figured recently that MATLAB's butter function was giving slightly different coefficients than those I calculate for a simple first order low pass filter, and I am not so sure why.

Here's how I calculate, starting from the following definition of a simple low passe filter and Bilinear Transform.

$$G(s)=\frac{1}{1+\frac{s}{\Omega_c}}$$ After bilinear transform, I get

$$ H(z) = G\bigg(\frac{2}{T}\frac{z-1}{z+1}\bigg) = \frac{1}{1+\frac{\frac{2}{T}\frac{z-1}{z+1}}{\Omega_c}} = \frac{\Omega_c}{\Omega_c+\frac{2}{T}\frac{z-1}{z+1}}$$ Then $$ H(z) = \frac{T\Omega_c(z+1)}{T\Omega_c(z+1)+2(z-1)} = \frac{T\Omega_c(z+1)}{T\Omega_cz+T\Omega_c+2z-2} = \frac{T\Omega_c(z+1)}{(T\Omega_c+2)z+(T\Omega_c-2)}$$ From there I define a substitution variable $k=\frac{T\Omega_c}{2}$ $$H(z) = \frac{2k(z+1)}{(2k+2)z+(2k-2)} = \frac{k(1+z^{-1})}{(k+1)+(k-1)z^{-1}}$$

Now, I compare my results with MATLAB's butter like this:

function fof( fc,fs)
    wc = 2*pi*fc;
    k = wc/(2*fs);
    b = [k k];
    a = [(k+1) (k-1)];
    b = b./a(1);    % Normalize coefficients for comparison only.
    a = a./a(1);

    [b2,a2] = butter(1,2*fc/fs);
    freqz(b,a,10*fs,fs);
    hold on
    freqz(b2,a2,10*fs,fs);
end

The bigger the cutoff frequency, the bigger the difference. It looks like MATLAB is indeed right as the freqz always shows a $-3\textrm{ dB}$ gain for the butter filter, and not mine.

Where am I wrong ?