Extract principal components
There is a nice trick to get the eigenvalues and eigenvectors of the covariance matrix without ever forming the matrix itself. This can be done using the singular value  decomposition (SVD), as described in this post from Stats.SE. Not only is this more numerically stable, but the results are automatically sorted.

A python version might look like this:

def components(X):
    _, vals, vecs = np.linalg.svd(X - X.mean(axis=0), full_matrices=False)
    return vals**2/(len(X)-1), vecs


A few things to note:


As described in the linked post above, the data matrix is typically defined with dimensions as columns, i.e. the transpose of your X.
The principle values and components are typically sorted from largest to smallest, i.e. the reverse of yours.
The function above does not assume that X has been pre-centered.


So to get results comparable to yours, you would need to do:

vals, vecs = components(X.T)
e_vals, e_vecs = vals[::-1], vecs[::-1]