diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx
index a5cfab6aa02..61993a04a27 100644
--- a/src/sage/matrix/matrix2.pyx
+++ b/src/sage/matrix/matrix2.pyx
@@ -1114,18 +1114,54 @@ cdef class Matrix(Matrix1):
             [-236774176922867   -3334450741532  470910201757143          1961587  230292926737068]
             [  82318322106118    1159275026338 -163719448527234          -681977  -80065012022313]
             [  53148766104440     748485096017 -105705345467375          -440318  -51693918051894]
+
+        TESTS:
+
+        Check for :issue:`40210`::
+
+            sage: A = vector(ZZ, [1, 2, 3]).column()
+            sage: B = vector(ZZ, [1, 1, 1]).column()
+            sage: A._solve_right_smith_form(B)
+            Traceback (most recent call last):
+            ...
+            ValueError: matrix equation has no solution
+
+        Random testing::
+
+            sage: n = randrange(1,100)
+            sage: m = randrange(1,100)
+            sage: A = matrix(ZZ, [[randrange(-10,11) for _ in range(m)] for _ in range(n)])
+            sage: y = A * vector(ZZ, [randrange(-100,101) for _ in range(m)])
+            sage: unsolvable = randrange(1 + (A.column_space() != ZZ^n))
+            sage: if unsolvable:
+            ....:     while y in A.column_space():
+            ....:         y += (ZZ^n).random_element()
+            sage: y = y.column()
+            sage: try:
+            ....:     x = A._solve_right_smith_form(y)
+            ....:     solved = True
+            ....: except ValueError:
+            ....:     solved = False
+            sage: solved == (not unsolvable)
+            True
+            sage: not solved or A * x == y
+            True
         """
         S,U,V = self.smith_form()
 
         n,m = self.dimensions()
         r = B.ncols()
 
+        UB = U * B
+        if UB[m:]:
+            raise ValueError("matrix equation has no solution")
+
         X_ = []
-        for d, v in zip(S.diagonal(), (U*B).rows()):
+        for d, v in zip(S.diagonal(), UB):
             if d:
                 X_.append(v / d)
             elif v:
-                raise ValueError("matrix equation has no solutions")
+                raise ValueError("matrix equation has no solution")
             else:
                 X_.append([0] * r)
 
@@ -1135,7 +1171,7 @@ cdef class Matrix(Matrix1):
         try:
             X_ = matrix(self.base_ring(), m, r, X_)
         except TypeError:
-            raise ValueError("matrix equation has no solutions")
+            raise ValueError("matrix equation has no solution")
 
         return V * X_
 
@@ -1181,26 +1217,65 @@ cdef class Matrix(Matrix1):
             [  968595469303  1461570161933   781069571508  1246248350502       -1629017]
             [ -235552378240  -355438713600  -189948023680  -303074680960         396160]
             [             0              0              0              0              0]
+
+        TESTS:
+
+        Check for :issue:`40210`::
+
+            sage: A = vector(ZZ, [1, 2, 3]).column()
+            sage: B = vector(ZZ, [1, 1, 1]).column()
+            sage: A._solve_right_hermite_form(B)
+            Traceback (most recent call last):
+            ...
+            ValueError: matrix equation has no solution
+
+        Random testing::
+
+            sage: n = randrange(1,100)
+            sage: m = randrange(1,100)
+            sage: A = matrix(ZZ, [[randrange(-10,11) for _ in range(m)] for _ in range(n)])
+            sage: y = A * vector(ZZ, [randrange(-100,101) for _ in range(m)])
+            sage: unsolvable = randrange(1 + (A.column_space() != ZZ^n))
+            sage: if unsolvable:
+            ....:     while y in A.column_space():
+            ....:         y += (ZZ^n).random_element()
+            sage: y = y.column()
+            sage: try:
+            ....:     x = A._solve_right_hermite_form(y)
+            ....:     solved = True
+            ....: except ValueError:
+            ....:     solved = False
+            sage: solved == (not unsolvable)
+            True
+            sage: not solved or A * x == y
+            True
         """
         H,U = self.transpose().hermite_form(transformation=True)
         H = H.transpose()
         U = U.transpose()
-#        assert self*U == H
 
         n,m = self.dimensions()
         r = B.ncols()
 
         from sage.matrix.constructor import matrix
         X_ = matrix(self.base_ring(), m, r)
-        for i in range(min(n,m)):
-            v = B[i,:]
-            v -= H[i,:i] * X_[:i]
-            d = H[i][i]
-            try:
-                X_[i] = v / d
-            except (ZeroDivisionError, TypeError) as e:
-                raise ValueError("matrix equation has no solution")
-#        assert H*X_ == B
+        j = 0  # current column
+        for i in range(n):
+            if j < m and H[i,j]:
+                # pivot for column j is in row i
+                v = B[i,:]
+                v -= H[i,:j] * X_[:j]
+                if v:
+                    try:
+                        X_[j] = v / H[i,j]
+                    except TypeError:
+                        raise ValueError("matrix equation has no solution")
+                j += 1
+            else:
+                # pivot for column j is below row i
+                assert not H[i,j:]
+                if H[i] * X_ != B[i]:
+                    raise ValueError("matrix equation has no solution")
 
         return U * X_