diff --git a/rrompy/utilities/poly_fitting/radial_basis/vander.py b/rrompy/utilities/poly_fitting/radial_basis/vander.py
index e798b23..90b3f2d 100644
--- a/rrompy/utilities/poly_fitting/radial_basis/vander.py
+++ b/rrompy/utilities/poly_fitting/radial_basis/vander.py
@@ -1,78 +1,78 @@
 # Copyright (C) 2018 by the RROMPy authors
 #
 # This file is part of RROMPy.
 #
 # RROMPy is free software: you can redistribute it and/or modify
 # it under the terms of the GNU Lesser General Public License as published by
 # the Free Software Foundation, either version 3 of the License, or
 # (at your option) any later version.
 #
 # RROMPy is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 # GNU Lesser General Public License for more details.
 #
 # You should have received a copy of the GNU Lesser General Public License
 # along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
 #
 
 import numpy as np
 from rrompy.utilities.poly_fitting.polynomial.vander import polyvander as pvP
 from rrompy.utilities.base.types import Np1D, Np2D, List, paramList
 from rrompy.parameter import checkParameterList
 from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
 from .kernel import radialGaussian, thinPlateSpline, multiQuadric
 
 __all__ = ['rbvander', 'polyvander']
 
 def rbvander(x:paramList, basis:str, reorder : List[int] = None,
              directionalWeights : Np1D = None) -> Np2D:
     """Compute radial-basis-Vandermonde matrix."""
     x, _ = checkParameterList(x)
     x_un, idx_un = x.unique(return_inverse = True)
     nx = len(x)
     if len(x_un) < nx:
         raise RROMPyException("Sample points must be distinct.")
     del x_un
     if directionalWeights is None:
         directionalWeights = np.ones(x.shape[1])
     RROMPyAssert(len(directionalWeights), x.shape[1],
                  "Number of directional weights")
     try:
         radialkernel = {"GAUSSIAN" : radialGaussian,
                         "THINPLATE" : thinPlateSpline,
                         "MULTIQUADRIC" : multiQuadric
                        }[basis.upper()]
     except:
         raise RROMPyException("Radial basis not recognized.")
     r2 = np.zeros((nx * (nx - 1) // 2 + 1), dtype = float)
     idxInv = np.zeros(nx ** 2, dtype = int)
     for j in range(nx):
         idx = j * (nx - 1) - j * (j + 1) // 2
         for i in range(j + 1, nx):
             r2[idx + i] = np.sum(np.abs((x[j] - x[i]) * directionalWeights
                                           ) ** 2.)
             idxInv[j * nx + i] = idx + i
             idxInv[i * nx + j] = idx + i
     Van = radialkernel(r2)[idxInv].reshape((nx, nx))
     if reorder is not None: Van = Van[reorder, :]
     return Van
     
 def polyvander(x:paramList, degs:List[int], basis:str,
                derIdxs : List[List[List[int]]] = None,
                reorder : List[int] = None, directionalWeights : Np1D = None,
                scl : Np1D = None) -> Np2D:
     """
     Compute radial-basis-inclusive Hermite-Vandermonde matrix with specified
         derivative directions.
     """
     if derIdxs is not None and np.sum(np.sum(derIdxs)) > 0:
         raise RROMPyException(("Cannot take derivatives of radial basis "
                                "function."))
     basisp, basisr = basis.split("_")
     VanR = rbvander(x, basisr, reorder = reorder,
                     directionalWeights = directionalWeights)
     VanP = pvP(x, degs, basisp, derIdxs = derIdxs, reorder = reorder,
                scl = scl)
-    return np.bmat([[VanR, VanP],
-                    [VanP.T.conj(), np.zeros(tuple([VanP.shape[1]] * 2))]])
+    return np.block([[VanR, VanP],
+                     [VanP.T.conj(), np.zeros(tuple([VanP.shape[1]] * 2))]])