docs: def block_diagonal_matrix(matrices, type=None):
      ur"""Build a block-diagonal matrix out of a given list of matrices.
  
      The type of matrix is chosen according to the input matrices.
  
      >>> import numpy as np
      >>> import sympy as sy
      >>> lis = [sy.ones(2), 2*sy.ones(3), 3*sy.ones(4)]
      >>> sy.pprint(block_diagonal_matrix(lis))
      ⎡1  1  0  0  0  0  0  0  0⎤
      ⎢                         ⎥
      ⎢1  1  0  0  0  0  0  0  0⎥
      ⎢                         ⎥
      ⎢0  0  2  2  2  0  0  0  0⎥
      ⎢                         ⎥
      ⎢0  0  2  2  2  0  0  0  0⎥
      ⎢                         ⎥
      ⎢0  0  2  2  2  0  0  0  0⎥
      ⎢                         ⎥
      ⎢0  0  0  0  0  3  3  3  3⎥
      ⎢                         ⎥
      ⎢0  0  0  0  0  3  3  3  3⎥
      ⎢                         ⎥
      ⎢0  0  0  0  0  3  3  3  3⎥
      ⎢                         ⎥
      ⎣0  0  0  0  0  3  3  3  3⎦
  
      >>> lis = [np.ones((2, 2)), 2*np.ones((3, 3)), 3*np.ones((4, 4))]
      >>> print(block_diagonal_matrix(lis))
      [[1. 1. 0. 0. 0. 0. 0. 0. 0.]
       [1. 1. 0. 0. 0. 0. 0. 0. 0.]
       [0. 0. 2. 2. 2. 0. 0. 0. 0.]
       [0. 0. 2. 2. 2. 0. 0. 0. 0.]
       [0. 0. 2. 2. 2. 0. 0. 0. 0.]
       [0. 0. 0. 0. 0. 3. 3. 3. 3.]
       [0. 0. 0. 0. 0. 3. 3. 3. 3.]
       [0. 0. 0. 0. 0. 3. 3. 3. 3.]
       [0. 0. 0. 0. 0. 3. 3. 3. 3.]]
  
      """
      if type is None:
          type = np.float64
      sizes = [Ai.shape[0] for Ai in matrices]
      size = sum(sizes)
      symbolic = hasattr(matrices[0][0], "subs")
      if symbolic:
          A = symzeros(size, size)
      else:
          A = np.zeros((size, size), type)
      ini = 0; fin = 0
      for i, sizei in enumerate(sizes):
          fin += sizei
          A[ini: fin, ini: fin] = matrices[i]
          ini += sizei
  
      return A