As of NumPy 1.13, there's numpy.block:

Block = numpy.block([[A11, A12], [A21, A22]])

For previous versions, there's bmat:

Block = numpy.bmat([[A11, A12], [A21, A22]])

numpy.bmat creates a matrix, rather than an array. This is usually a bad thing. You can call asarray on the result if you want an array, or use the A attribute:

Block = numpy.bmat([[A11, A12], [A21, A22]]).A

bmat also does some messing around with stack frames to let you do this:

Block = numpy.bmat('A11,A12; A21,A22')

Here are some examples of a different (loop free) way to work with blocks:

import numpy as np
from numpy.lib.stride_tricks import as_strided as ast

A= np.arange(36).reshape(6, 6)
print A
#[[ 0  1  2  3  4  5]
# [ 6  7  8  9 10 11]
# ...
# [30 31 32 33 34 35]]

# 2x2 block view
B= ast(A, shape= (3, 3, 2, 2), strides= (48, 8, 24, 4))
print B[1, 1]
#[[14 15]
# [20 21]]

# for preserving original shape
B[:, :]= np.dot(B[:, :], np.array([[0, 1], [1, 0]]))
print A
#[[ 1  0  3  2  5  4]
# [ 7  6  9  8 11 10]
# ...
# [31 30 33 32 35 34]]
print B[1, 1]
#[[15 14]
# [21 20]]

# for reducing shape, processing in 3D is enough
C= B.reshape(3, 3, -1)
print C.sum(-1)
#[[ 14  22  30]
# [ 62  70  78]
# [110 118 126]]

So just trying to simply copy the matlab functionality to numpy is not all ways the best way to proceed. Sometimes a 'off the hat' thinking is needed.

Caveat:
In general, implementations based on stride tricks may (but does not necessary need to) suffer some performance penalties. So be prepared to all ways measure your performance. In any case it's wise to first check if the needed functionality (or similar enough, in order to easily adapt for) has all ready been implemented in numpy or scipy.

Update:
Please note that there is no real magic involved here with the strides, so I'll provide a simple function to get a block_view of any suitable 2D numpy-array. So here we go:

from numpy.lib.stride_tricks import as_strided as ast

def block_view(A, block= (3, 3)):
    """Provide a 2D block view to 2D array. No error checking made.
    Therefore meaningful (as implemented) only for blocks strictly
    compatible with the shape of A."""
    # simple shape and strides computations may seem at first strange
    # unless one is able to recognize the 'tuple additions' involved ;-)
    shape= (A.shape[0]/ block[0], A.shape[1]/ block[1])+ block
    strides= (block[0]* A.strides[0], block[1]* A.strides[1])+ A.strides
    return ast(A, shape= shape, strides= strides)

if __name__ == '__main__':
    from numpy import arange
    A= arange(144).reshape(12, 12)
    print block_view(A)[0, 0]
    #[[ 0  1  2]
    # [12 13 14]
    # [24 25 26]]
    print block_view(A, (2, 6))[0, 0]
    #[[ 0  1  2  3  4  5]
    # [12 13 14 15 16 17]]
    print block_view(A, (3, 12))[0, 0]
    #[[ 0  1  2  3  4  5  6  7  8  9 10 11]
    # [12 13 14 15 16 17 18 19 20 21 22 23]
    # [24 25 26 27 28 29 30 31 32 33 34 35]]

The code proposed by the OP can indeed made be more efficient, mainly by noting the fact that to form the sequence $A^i B$, with $i=0\,\dots,N$ you do not have to compute $A^i$ at each step, but you can exploit the fact that $A^i B = A\,(A^{i-1}B)$, reusing the result of the previous step.

My proposed implementation is

import numpy as np

N = 3
M = 1

A = np.random.random((N, N))
B = np.random.random((N, M))

X = np.zeros((N ** 2, N * M))
rsl = slice(0, N)
X[rsl, :M] = B
for i in range(1, N):
    rsl_p, rsl = rsl, slice(i * N, (i + 1) * N)
    X[rsl, :M] = A @ X[rsl_p, :M]
    X[rsl, M : (i + 1) * M] = X[rsl_p, : i * M]

This question is more related to programming that computational science, so I add some general remarks, useful for a novice programmer.

1. never blindly translate mathematical expressions into code: think algorithmically.

   X = np.linalg.matrix_power(A,i)@B 
   

   is much slower than

   X=B; 
   for j in range(i): X = A*X
   

   if $M \ll N$

   (exercise: compute operation count for both code fragments).

2. whenever possible reuse previous calculations.

3. pre-allocate output matrices, and do not store intermediate results in temp. matrices

mats_row=np.array([[A1,A2,A3,...,A_n]]) #Create array of matrices with shape (1,M,N,N)
mats_column=np.transpose(mats_row,(1,0,2,3)) #Make a copy with shape (M,1,N,N)
block=.5*(mats_row+mats_column) #Add, broadcasting to a (M,M,N,N) array

This works by utilizing Numpy's broadcasting capabilities. The basic idea is similar to if you added a matrix where each row was $[A_1,A_2,..A_n]$ to a matrix where each column was $[A_1,A_2,..A_n]$, but broadcasting allows you to do this without ever explicitly forming these intermediates matrices.

IMO, np.eye already has everything you need, as you can define number of rows and columns separately.
So your function should simply look like

def fct(k):
    return np.eye(k**2-k, k**2)

Easy way to represent your matrix in blocks (keeping intermediate steps, but you don't have to):

block_m = 8  # rows in block
block_n = 8  # cols in block
A = np.arange(256*256).reshape(256,256)  # sample data
B = A.reshape((-1, block_m, A.shape[1]//block_n, block_n))  
C = B.transpose((0,2,1,3))

This yields (in this case) a 32-by-32-by-8-by-8 array where the first two axes refer to the block matrix. i.e. C[0,0] is the 8-by-8 sub-array in the top left of your original array. In most cases, I believe C should be a view on the original array; changing C[0,0,0,0], for instance, changes A[0,0].

You can loop over the first two axes or apply vectorized functions as needed. As an example, here's me averaging the 8-by-8 subarrays:

block_means = np.sum(C, axis=(2,3)) / 64