cov(X,Y) is equivalent to cov([x(:) y(:)]). But [x(:) y(:)] is 20000 by 2 for you, and cov() treats rows as observations and columns as dimensions, so you get a 2 by 2 covariance matrix.

I would just implement it myself by the definition:

bsxfun(@minus,x,mean(x))'*bsxfun(@minus,y,mean(y))/(size(x,1)-1)

If you have an older version of matlab that doesn't support bsxfun(), just use repmat().

Note that numpy.cov() considers its input data matrix to have observations in each column, and variables in each row, so to get numpy.cov() to return what other packages do, you have to pass the transpose of the data matrix to numpy.cov().

The Python code that you linked can be used to simulate what other packages do, but it contains some errors: N should be the number of rows, not columns, and you have to perform the matrix multiplication in the other order:

import numpy as np

def cov(X0):
    print "\n==\nMatrix:"
    print X0
    X = X0 - X0.mean(axis=0)
    N = X.shape[0]                # !!!
    fact = float(N - 1)
    print "Covariance:"
    print np.dot(X.T, X) / fact   # !!!

X0 = np.vstack(([1, 2], [3, 4]))
cov(X0)
cov(X0.T)

X0 = np.vstack(([1, 2], [3, 4], [22, 44]))
cov(X0)
cov(X0.T)

With these fixes, the covariance behaves as expected:

==
Matrix:
[[1 2]
 [3 4]]
Covariance:
[[ 2.  2.]
 [ 2.  2.]]

==
Matrix:
[[1 3]
 [2 4]]
Covariance:
[[ 0.5  0.5]
 [ 0.5  0.5]]

==
Matrix:
[[ 1  2]
 [ 3  4]
 [22 44]]
Covariance:
[[ 134.33333333  274.33333333]
 [ 274.33333333  561.33333333]]

==
Matrix:
[[ 1  3 22]
 [ 2  4 44]]
Covariance:
[[   0.5    0.5   11. ]
 [   0.5    0.5   11. ]
 [  11.    11.   242. ]]

As for the numpy.linalg.svd() code, you need to center the data matrix by subtracting off the variable means, and the multiplication involving the V matrix must be performed in the other order. With these changes you will replicate everybody else's behavior:

import numpy as np

def cov(X0):
    print "\n==\nMatrix:"
    print X0
    X = X0 - X0.mean(axis=0)
    U, s, V = np.linalg.svd(X, full_matrices = 0)
    D = np.dot(np.dot(V.T,np.diag(s**2)),V)
    Dadjust = D / (X0.shape[0] - 1)
    print "Covariance:"
    print (Dadjust)