diff --git a/pyproject.toml b/pyproject.toml
index 470323e..2759790 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,6 +1,6 @@
 [project]
 name = "calfem-python"
-version = "3.6.12"
+version = "3.6.13"
 description = "CALFEM for Python"
 authors = [
     {name = "Jonas Lindemann", email = "jonas.lindemann@lunarc.lu.se"},
diff --git a/src/calfem/core.py b/src/calfem/core.py
index 5e0b946..52b627b 100644
--- a/src/calfem/core.py
+++ b/src/calfem/core.py
@@ -5473,7 +5473,7 @@ def solveq(K, f, bcPrescr=None, bcVal=None):
     """
 
     if bcPrescr is None:
-        return np.asmatrix(np.linalg.solve(K, f))
+        return np.array(np.linalg.solve(K, f))
 
     nDofs = K.shape[0]
     nPdofs = bcPrescr.shape[0]
@@ -5487,17 +5487,20 @@ def solveq(K, f, bcPrescr=None, bcVal=None):
     bc[np.ix_(bcPrescr-1)] = False
     bcDofs = bcDofs[bc]
 
-    fsys = f[bcDofs]-K[np.ix_((bcDofs), (bcPrescr-1))] * \
-        np.asmatrix(bcVal).reshape(nPdofs, 1)
+    # Ensure bcVal is a column vector
+    bcVal_col = np.array(bcVal).reshape(-1, 1)
+
+    # Compute fsys with correct broadcasting
+    fsys = f[bcDofs] - K[np.ix_((bcDofs), (bcPrescr-1))] @ bcVal_col
     asys = np.linalg.solve(K[np.ix_((bcDofs), (bcDofs))], fsys)
 
     a = np.zeros([nDofs, 1])
-    a[np.ix_(bcPrescr-1)] = np.asmatrix(bcVal).reshape(nPdofs, 1)
-    a[np.ix_(bcDofs)] = asys
+    a[np.ix_(bcPrescr-1)] = bcVal_col
+    a[np.ix_(bcDofs)] = asys.reshape(-1, 1)
 
-    Q = K*np.asmatrix(a)-f
+    Q = K @ a - f
 
-    return (np.asmatrix(a), Q)
+    return (a, Q)
 
 
 def spsolveq(K, f, bcPrescr, bcVal=None):
@@ -5536,7 +5539,7 @@ def spsolveq(K, f, bcPrescr, bcVal=None):
     bc[np.ix_(bcPrescr-1)] = False
     bcDofs = bcDofs[bc]
 
-    bcVal_m = np.asmatrix(bcVal).reshape(nPdofs, 1)
+    bcVal_m = np.array(bcVal).reshape(-1, 1)
 
     info("Preparing system matrix...")
 
@@ -5550,7 +5553,7 @@ def spsolveq(K, f, bcPrescr, bcVal=None):
     #Kt = Kt1[:,bcPrescr]
     Kt = K[np.ix_((bcDofs), (bcPrescr-1))]
     info("step 3... fsys")
-    fsys = f[bcDofs]-Kt*bcVal_m
+    fsys = f[bcDofs] - Kt @ bcVal_m
     info("step 4... Ksys")
     Ksys1 = Kcsr[bcDofs]
     Ksys = Ksys1[:, bcDofs]
@@ -5563,12 +5566,12 @@ def spsolveq(K, f, bcPrescr, bcVal=None):
     info("Reconstructing full a...")
     a = np.zeros([nDofs, 1])
     a[np.ix_(bcPrescr-1)] = bcVal_m
-    a[np.ix_(bcDofs)] = np.asmatrix(asys).transpose()
+    # asys is 1D, ensure shape (n,1) for assignment
+    a[np.ix_(bcDofs)] = np.array(asys).reshape(-1, 1)
 
-    a_m = np.asmatrix(a)
-    Q = K*a_m-f
+    Q = K @ a - f
     info("done...")
-    return (a_m, Q)
+    return (a, Q)
 
 
 def eigen(K: ArrayLike, M: ArrayLike, b: Optional[ArrayLike] = None) -> Tuple[NDArray[np.floating], NDArray[np.floating]]:
diff --git a/src/calfem/matrix_compat.py b/src/calfem/matrix_compat.py
index 7c9b6a3..901b702 100644
--- a/src/calfem/matrix_compat.py
+++ b/src/calfem/matrix_compat.py
@@ -299,6 +299,15 @@ def flatten(self):
         Return a flattened copy of the array
         """
         return self.array.flatten()
+    
+def asmatrix(data):
+    """
+    Compatibility replacement for np.asmatrix.
+    Converts input to MatrixCompat if not already.
+    """
+    if isinstance(data, MatrixCompat):
+        return data
+    return MatrixCompat(data)
 
 # Now we patch numpy to use our compatibility layer, but only for NumPy >= 2.0
 # where np.matrix has been removed
@@ -314,9 +323,12 @@ def _get_numpy_version():
     np_matrix_original = getattr(np, 'matrix', None)
     np.matrix = MatrixCompat
     np.mat = MatrixCompat  # Alias for np.matrix
+    np_asmatrix_original = getattr(np, 'asmatrix', None)
+    np.asmatrix = asmatrix
 else:
     # For NumPy < 2.0, keep the original matrix
     np_matrix_original = np.matrix
+    np_asmatrix_original = np.asmatrix
 
 # Optional cleanup function to restore original np.matrix if needed
 def restore_numpy_matrix():
@@ -324,4 +336,5 @@ def restore_numpy_matrix():
     Restore the original np.matrix function
     """
     np.matrix = np_matrix_original
-    np.mat = np_matrix_original  # Restore alias as well
\ No newline at end of file
+    np.mat = np_matrix_original  # Restore alias as well
+    np.asmatrix = np_asmatrix_original
\ No newline at end of file