diff --git a/smErrors/dimension_error.go b/smErrors/dimension_error.go
index 419a16e..fc69d8f 100644
--- a/smErrors/dimension_error.go
+++ b/smErrors/dimension_error.go
@@ -125,3 +125,71 @@ func CheckDimensionsInMultiplication(left, right MatrixLike) error {
 	// If dimensions match, then return nothing.
 	return nil
 }
+
+/*
+CheckDimensionsInHStack
+Description:
+
+	This function checks that the dimensions of the left and right expressions
+	are compatible for horizontal stacking.
+	We allow:
+	- Stacking if the number of rows match
+*/
+func CheckDimensionsInHStack(sliceToStack ...MatrixLike) error {
+	// Check that the size of columns in left and right agree
+	var nRowsInSlice []int
+	for _, slice := range sliceToStack {
+		nRowsInSlice = append(nRowsInSlice, slice.Dims()[0])
+	}
+
+	// Check that the number of rows in each slice is the same
+	for ii := 1; ii < len(nRowsInSlice); ii++ {
+		// If the number of rows in the slice is not the same as the previous slice,
+		// then return an error
+		dimsAreMatched := nRowsInSlice[ii] == nRowsInSlice[ii-1]
+		if !dimsAreMatched {
+			return DimensionError{
+				Operation: "HStack",
+				Arg1:      sliceToStack[ii-1],
+				Arg2:      sliceToStack[ii],
+			}
+		}
+	}
+
+	// If dimensions match, then return nothing.
+	return nil
+}
+
+/*
+CheckDimensionsInVStack
+Description:
+
+	This function checks that the dimensions of the left and right expressions
+	are compatible for vertical stacking.
+	We allow:
+	- Stacking if the number of columns match
+*/
+func CheckDimensionsInVStack(sliceToStack ...MatrixLike) error {
+	// Check that the size of columns in left and right agree
+	var nColsInSlice []int
+	for _, slice := range sliceToStack {
+		nColsInSlice = append(nColsInSlice, slice.Dims()[1])
+	}
+
+	// Check that the number of rows in each slice is the same
+	for ii := 1; ii < len(nColsInSlice); ii++ {
+		// If the number of rows in the slice is not the same as the previous slice,
+		// then return an error
+		dimsAreMatched := nColsInSlice[ii] == nColsInSlice[ii-1]
+		if !dimsAreMatched {
+			return DimensionError{
+				Operation: "VStack",
+				Arg1:      sliceToStack[ii-1],
+				Arg2:      sliceToStack[ii],
+			}
+		}
+	}
+
+	// If dimensions match, then return nothing.
+	return nil
+}
diff --git a/symbolic/constant.go b/symbolic/constant.go
index babd937..6d08e01 100644
--- a/symbolic/constant.go
+++ b/symbolic/constant.go
@@ -70,7 +70,7 @@ func (c K) LinearCoeff(wrt ...[]Variable) mat.VecDense {
 		// If the user didn't provide any variables, then panic!
 		// We cannot construct zero length vectors in gonum
 		panic(
-			smErrors.EmptyLinearCoeffsError{c},
+			smErrors.EmptyLinearCoeffsError{Expression: c},
 		)
 	}
 
diff --git a/symbolic/constant_vector.go b/symbolic/constant_vector.go
index f19b048..638aedf 100644
--- a/symbolic/constant_vector.go
+++ b/symbolic/constant_vector.go
@@ -75,7 +75,7 @@ func (kv KVector) AtVec(idx int) ScalarExpression {
 	// Input Processing
 
 	// Check to see whether or not the index is valid.
-	err := smErrors.CheckIndexOnMatrix(idx, 0, mc)
+	err := smErrors.CheckIndexOnMatrix(idx, 0, kv)
 	if err != nil {
 		panic(err)
 	}
diff --git a/symbolic/expression.go b/symbolic/expression.go
index e39823d..52407ed 100644
--- a/symbolic/expression.go
+++ b/symbolic/expression.go
@@ -4,7 +4,6 @@ import (
 	"fmt"
 
 	"github.com/MatProGo-dev/SymbolicMath.go/smErrors"
-	"github.com/MatProGo-dev/SymbolicMath.go/symbolic"
 )
 
 /*
@@ -70,7 +69,7 @@ type Expression interface {
 	Power(exponent int) Expression
 
 	// At returns the value at the given row and column index
-	At(rowIndex, colIndex, nCols int) symbolic.ScalarExpression
+	At(ii, jj int) ScalarExpression
 }
 
 /*
@@ -215,18 +214,118 @@ func HStack(eIn ...Expression) Expression {
 	}
 
 	// Create the resulting Matrix's shape
-	var result [][]symbolic.ScalarExpression
+	var result [][]ScalarExpression
 	for rowIndex := 0; rowIndex < eIn[0].Dims()[0]; rowIndex++ {
-		var tempRow []symbolic.ScalarExpression
-		for stackIndex_ii := 0; stackIndex_ii < len(eIn); stackIndex_ii++ {
-			nCols_ii := nCols[stackIndex_ii]
+		var tempRow []ScalarExpression
+		for stackIndexII := 0; stackIndexII < len(eIn); stackIndexII++ {
+			nCols_ii := nCols[stackIndexII]
 			// Add all of the columns from the current expression to the row
 			for colIndex := 0; colIndex < nCols_ii; colIndex++ {
-				tempRow = append(tempRow, eIn[stackIndex_ii].At(rowIndex, colIndex, nCols_ii))
+				tempRow = append(tempRow, eIn[stackIndexII].At(rowIndex, colIndex))
 			}
 		}
+		// Add the row to the result
+		result = append(result, tempRow)
 	}
 
 	// Return the simplified form of the expression
-	return symbolic.ConcretizeMatrixExpression(result)
+	return ConcretizeExpression(result)
+}
+
+/*
+VStack
+Description:
+
+	Stacks the input expressions vertically.
+*/
+func VStack(eIn ...Expression) Expression {
+	// Input Checking
+
+	// Panic if there are 0 expressions in the input
+	if len(eIn) == 0 {
+		panic(
+			fmt.Errorf("VStack: There must be at least one expression in the input; received 0"),
+		)
+	}
+
+	// Check that all the expressions have the same number of columns
+	var mlSlice []smErrors.MatrixLike // First convert expression slice to matrix like slice
+	for _, e := range eIn {
+		mlSlice = append(mlSlice, e)
+	}
+
+	err := smErrors.CheckDimensionsInVStack(mlSlice...)
+	if err != nil {
+		panic(err)
+	}
+
+	// Setup
+
+	// Create the resulting Matrix's shape
+	var result [][]ScalarExpression
+	for stackIndexII := 0; stackIndexII < len(eIn); stackIndexII++ {
+		// Each row will be made from the rows of the current expression
+		eII := eIn[stackIndexII]
+		for rowIndex := 0; rowIndex < eII.Dims()[0]; rowIndex++ {
+			var tempRow []ScalarExpression
+			for colIndex := 0; colIndex < eII.Dims()[1]; colIndex++ {
+				tempRow = append(tempRow, eII.At(rowIndex, colIndex))
+			}
+			// Add the row to the result
+			result = append(result, tempRow)
+		}
+	}
+
+	// Return the simplified form of the expression
+	return ConcretizeExpression(result)
+}
+
+/*
+ConcretizeExpression
+Description:
+
+	Converts the input expression to a valid type that implements "Expression".
+*/
+func ConcretizeExpression(e interface{}) Expression {
+	// Input Processing
+
+	// Convert
+	var (
+		concrete Expression
+	)
+	switch e.(type) {
+	case []ScalarExpression:
+		concreteVectorE := ConcretizeVectorExpression(e.([]ScalarExpression))
+		// If vector expression is a scalar (i.e., has 1 row), return the scalar expression
+		if concreteVectorE.Dims()[0] == 1 {
+			concrete = concreteVectorE.At(0, 0)
+		} else {
+			concrete = concreteVectorE
+		}
+
+	case [][]ScalarExpression:
+		concreteMatrixE := ConcretizeMatrixExpression(e.([][]ScalarExpression))
+		// If matrix expression is a scalar (i.e., has 1 row and 1 column), return the scalar expression
+		switch {
+		case concreteMatrixE.Dims()[0] == 1 && concreteMatrixE.Dims()[1] == 1: // If the matrix is a scalar
+			concrete = concreteMatrixE.At(0, 0)
+		case concreteMatrixE.Dims()[1] == 1: // If the matrix is a column vector
+			interm := make([]ScalarExpression, concreteMatrixE.Dims()[0])
+			for ii := 0; ii < concreteMatrixE.Dims()[0]; ii++ {
+				interm[ii] = concreteMatrixE.At(ii, 0)
+			}
+			concrete = ConcretizeVectorExpression(interm)
+		default:
+			concrete = concreteMatrixE
+		}
+	default:
+		panic(
+			smErrors.UnsupportedInputError{
+				FunctionName: "ConcretizeExpression",
+				Input:        e,
+			},
+		)
+	}
+
+	return concrete
 }
diff --git a/symbolic/matrix_expression.go b/symbolic/matrix_expression.go
index c6ead3a..ef74937 100644
--- a/symbolic/matrix_expression.go
+++ b/symbolic/matrix_expression.go
@@ -2,6 +2,7 @@ package symbolic
 
 import (
 	"fmt"
+
 	"github.com/MatProGo-dev/SymbolicMath.go/smErrors"
 	"gonum.org/v1/gonum/mat"
 )
@@ -353,7 +354,7 @@ func ConcretizeMatrixExpression(sliceIn [][]ScalarExpression) MatrixExpression {
 		return out
 
 	case isAllVariables:
-		// Convert to a variable vector
+		// Convert to a variable matrix
 		var out VariableMatrix
 		for _, row_ii := range sliceIn {
 			var tempRow []Variable
diff --git a/symbolic/monomial.go b/symbolic/monomial.go
index 70fb318..880b334 100644
--- a/symbolic/monomial.go
+++ b/symbolic/monomial.go
@@ -398,7 +398,7 @@ func (m Monomial) LinearCoeff(wrt ...[]Variable) mat.VecDense {
 	}
 
 	if len(wrtVars) == 0 {
-		panic(smErrors.CanNotGetLinearCoeffOfConstantError{m})
+		panic(smErrors.CanNotGetLinearCoeffOfConstantError{Expression: m})
 	}
 
 	// Algorithm
diff --git a/symbolic/monomial_matrix.go b/symbolic/monomial_matrix.go
index eb68549..dbe8f78 100644
--- a/symbolic/monomial_matrix.go
+++ b/symbolic/monomial_matrix.go
@@ -2,6 +2,7 @@ package symbolic
 
 import (
 	"fmt"
+
 	"github.com/MatProGo-dev/SymbolicMath.go/smErrors"
 	"gonum.org/v1/gonum/mat"
 )
@@ -35,7 +36,7 @@ Description:
 func (mm MonomialMatrix) Check() error {
 	// Check that the matrix has at least one row
 	if len(mm) == 0 {
-		return smErrors.EmptyMatrixError{mm}
+		return smErrors.EmptyMatrixError{Expression: mm}
 	}
 
 	// Check that the number of columns is the same in each row
diff --git a/symbolic/monomial_vector.go b/symbolic/monomial_vector.go
index af819ba..868f502 100644
--- a/symbolic/monomial_vector.go
+++ b/symbolic/monomial_vector.go
@@ -31,7 +31,7 @@ Description:
 func (mv MonomialVector) Check() error {
 	// Check that the polynomial has at least one monomial
 	if len(mv) == 0 {
-		return smErrors.EmptyVectorError{mv}
+		return smErrors.EmptyVectorError{Expression: mv}
 	}
 
 	// Check that each of the monomials are well formed
diff --git a/symbolic/polynomial.go b/symbolic/polynomial.go
index 0577f0d..a42bb01 100644
--- a/symbolic/polynomial.go
+++ b/symbolic/polynomial.go
@@ -702,7 +702,7 @@ func (p Polynomial) LinearCoeff(wrt ...[]Variable) mat.VecDense {
 
 	// If there are no variables in the slice, then return a vector of length 1 containing zero.
 	if len(wrtVars) == 0 {
-		panic(smErrors.CanNotGetLinearCoeffOfConstantError{p})
+		panic(smErrors.CanNotGetLinearCoeffOfConstantError{Expression: p})
 	}
 
 	// Algorithm
diff --git a/symbolic/polynomial_matrix.go b/symbolic/polynomial_matrix.go
index f22ed50..bfa3ea4 100644
--- a/symbolic/polynomial_matrix.go
+++ b/symbolic/polynomial_matrix.go
@@ -2,6 +2,7 @@ package symbolic
 
 import (
 	"fmt"
+
 	"github.com/MatProGo-dev/SymbolicMath.go/smErrors"
 	"gonum.org/v1/gonum/mat"
 )
@@ -35,7 +36,7 @@ Description:
 func (pm PolynomialMatrix) Check() error {
 	// Check that the matrix has at least one row
 	if len(pm) == 0 {
-		return smErrors.EmptyMatrixError{pm}
+		return smErrors.EmptyMatrixError{Expression: pm}
 	}
 
 	// Check that the number of columns is the same in each row
diff --git a/symbolic/polynomial_vector.go b/symbolic/polynomial_vector.go
index c2d0c20..67b8eab 100644
--- a/symbolic/polynomial_vector.go
+++ b/symbolic/polynomial_vector.go
@@ -181,7 +181,7 @@ func (pv PolynomialVector) LinearCoeff(vSlices ...[]Variable) mat.Dense {
 
 	if len(varSlice) == 0 {
 		panic(
-			smErrors.CanNotGetLinearCoeffOfConstantError{pv},
+			smErrors.CanNotGetLinearCoeffOfConstantError{Expression: pv},
 		)
 	}
 
diff --git a/symbolic/variable.go b/symbolic/variable.go
index e7e1c16..20bd91a 100644
--- a/symbolic/variable.go
+++ b/symbolic/variable.go
@@ -60,7 +60,7 @@ func (v Variable) LinearCoeff(wrt ...[]Variable) mat.VecDense {
 	}
 
 	if len(wrtVars) == 0 {
-		panic(smErrors.CanNotGetLinearCoeffOfConstantError{v})
+		panic(smErrors.CanNotGetLinearCoeffOfConstantError{Expression: v})
 	}
 
 	// Constants
diff --git a/symbolic/variable_matrix.go b/symbolic/variable_matrix.go
index 37e78ef..0e05f6d 100644
--- a/symbolic/variable_matrix.go
+++ b/symbolic/variable_matrix.go
@@ -91,9 +91,7 @@ func (vm VariableMatrix) Variables() []Variable {
 	// Algorithm
 	var variables []Variable
 	for _, vmRow := range vm {
-		for _, v := range vmRow {
-			variables = append(variables, v)
-		}
+		variables = append(variables, vmRow...) // Unrolls all of vmRow and appends it to variables
 	}
 
 	return UniqueVars(variables)
@@ -629,7 +627,7 @@ func NewVariableMatrixClassic(nRows, nCols int, envs ...*Environment) VariableMa
 		env = envs[0]
 	default:
 		panic(
-			fmt.Errorf("Too many inputs provided to NewVariableMatrix() method."),
+			fmt.Errorf("Too many inputs provided to NewVariableMatrix() method"),
 		)
 	}
 
diff --git a/testing/symbolic/expression_test.go b/testing/symbolic/expression_test.go
index 10de405..8d40b6b 100644
--- a/testing/symbolic/expression_test.go
+++ b/testing/symbolic/expression_test.go
@@ -1,8 +1,10 @@
 package symbolic_test
 
 import (
-	"github.com/MatProGo-dev/SymbolicMath.go/symbolic"
 	"testing"
+
+	"github.com/MatProGo-dev/SymbolicMath.go/smErrors"
+	"github.com/MatProGo-dev/SymbolicMath.go/symbolic"
 )
 
 /*
@@ -76,3 +78,389 @@ func TestExpression_ToExpression1(t *testing.T) {
 	}()
 	symbolic.ToExpression("x")
 }
+
+/*
+TestExpression_HStack1
+Description:
+
+	Tests the HStack function for two variables.
+	The result of the stacking should be a variable matrix with one row and two columns.
+*/
+func TestExpression_HStack1(t *testing.T) {
+	// Constants
+	x := symbolic.NewVariable()
+	y := symbolic.NewVariable()
+
+	// Test
+	result := symbolic.HStack(x, y)
+	if result.Dims()[0] != 1 || result.Dims()[1] != 2 {
+		t.Errorf(
+			"Expected the result to be a 1x2 matrix; received %v",
+			result.Dims(),
+		)
+	}
+
+	// Verify that the result is a variable matrix
+	if _, ok := result.(symbolic.VariableMatrix); !ok {
+		t.Errorf(
+			"Expected the result to be a VariableMatrix; received %T",
+			result,
+		)
+	}
+}
+
+/*
+TestExpression_HStack2
+Description:
+
+	Tests the HStack function for 4 scalar expressions. 3 of the expressions are
+	constants and the last one is a variable.
+	The result should be a monomial matrix with one row and 4 columns.
+*/
+func TestExpression_HStack2(t *testing.T) {
+	// Constants
+	x := symbolic.NewVariable()
+	c1 := symbolic.K(1.0)
+	c2 := symbolic.K(2.0)
+	c3 := symbolic.K(3.0)
+
+	// Test
+	result := symbolic.HStack(c1, c2, c3, x)
+	if result.Dims()[0] != 1 || result.Dims()[1] != 4 {
+		t.Errorf(
+			"Expected the result to be a 1x4 matrix; received %v",
+			result.Dims(),
+		)
+	}
+
+	// Verify that the result is a monomial matrix
+	if _, ok := result.(symbolic.MonomialMatrix); !ok {
+		t.Errorf(
+			"Expected the result to be a MonomialMatrix; received %T",
+			result,
+		)
+	}
+}
+
+/*
+TestExpression_HStack3
+Description:
+
+	Tests the HStack function for 2 vector expressions.
+	Each vector has 11 elements. One is a constant vector and the other is a variable vector.
+	The result should be a monomial matrix with 11 rows and 2 columns.
+*/
+func TestExpression_HStack3(t *testing.T) {
+	// Constants
+	kv1 := symbolic.VecDenseToKVector(symbolic.OnesVector(11))
+	vv2 := symbolic.NewVariableVector(11)
+
+	// Test
+	result := symbolic.HStack(kv1, vv2)
+	if result.Dims()[0] != 11 || result.Dims()[1] != 2 {
+		t.Errorf(
+			"Expected the result to be an 11x2 matrix; received %v",
+			result.Dims(),
+		)
+	}
+
+	// Verify that the result is a monomial matrix
+	if _, ok := result.(symbolic.MonomialMatrix); !ok {
+		t.Errorf(
+			"Expected the result to be a MonomialMatrix; received %T",
+			result,
+		)
+	}
+}
+
+/*
+TestExpression_HStack4
+Description:
+
+	Tests the HStack function for a matrix and a vector expression.
+	The matrix is a constant matrix and the vector is a variable vector.
+	The matrix is of shape 3x2 and the vector is of length 2.
+	The HStack function should panic because the dimensions do not match.
+*/
+func TestExpression_HStack4(t *testing.T) {
+	// Constants
+	km1 := symbolic.DenseToKMatrix(symbolic.OnesMatrix(3, 2))
+	vv2 := symbolic.NewVariableVector(2)
+
+	// Test
+	defer func() {
+		err := recover().(error)
+		if err == nil {
+			t.Errorf("The HStack function should panic when the dimensions do not match")
+		}
+
+		// Collect the expected error which should be a dimension error and
+		// compare it with the recovered error
+		expectedError := smErrors.DimensionError{
+			Operation: "HStack",
+			Arg1:      km1,
+			Arg2:      vv2,
+		}
+		if err.Error() != expectedError.Error() {
+			t.Errorf(
+				"Expected the error to be %v; received %v",
+				expectedError,
+				err,
+			)
+		}
+
+	}()
+	symbolic.HStack(km1, vv2)
+
+	// The function should panic before this point
+	t.Errorf("The HStack function should panic before this point")
+
+}
+
+/*
+TestExpression_HStack5
+Description:
+
+	Tests the HStack function for a matrix and a vector expression.
+	The matrix is a constant matrix and the vector is a variable vector.
+	The matrix is of shape 3x2 and the vector is of length 3.
+	The HStack function should not panic because the dimensions match.
+	It should produce a monomial matrix with 3 rows and 3 columns.
+*/
+func TestExpression_HStack5(t *testing.T) {
+	// Constants
+	km1 := symbolic.DenseToKMatrix(symbolic.OnesMatrix(3, 2))
+	vv3 := symbolic.NewVariableVector(3)
+
+	// Test
+	result := symbolic.HStack(km1, vv3)
+	if result.Dims()[0] != 3 || result.Dims()[1] != 3 {
+		t.Errorf(
+			"Expected the result to be a 3x3 matrix; received %v",
+			result.Dims(),
+		)
+	}
+
+	// Verify that the result is a monomial matrix
+	if _, ok := result.(symbolic.MonomialMatrix); !ok {
+		t.Errorf(
+			"Expected the result to be a MonomialMatrix; received %T",
+			result,
+		)
+	}
+}
+
+/*
+TestExpression_VStack1
+Description:
+
+	Tests the VStack function for two variables.
+	The result of the stacking should be a variable vector with two elements.
+*/
+func TestExpression_VStack1(t *testing.T) {
+	// Constants
+	x := symbolic.NewVariable()
+	y := symbolic.NewVariable()
+
+	// Test
+	result := symbolic.VStack(x, y)
+	if (result.Dims()[0] != 2) || (result.Dims()[1] != 1) {
+		t.Errorf(
+			"Expected the result to be a vector with 2 elements; received object with shape %v",
+			result.Dims(),
+		)
+	}
+
+	// Verify that the result is a variable vector
+	if _, ok := result.(symbolic.VariableVector); !ok {
+		t.Errorf(
+			"Expected the result to be a VariableVector; received %T",
+			result,
+		)
+	}
+}
+
+/*
+TestExpression_VStack2
+Description:
+
+	Tests the VStack function for 4 scalar expressions. 3 of the expressions are
+	constants and the last one is a variable.
+	The result should be a monomial vector with 4 elements.
+*/
+func TestExpression_VStack2(t *testing.T) {
+	// Constants
+	x := symbolic.NewVariable()
+	c1 := symbolic.K(1.0)
+	c2 := symbolic.K(2.0)
+	c3 := symbolic.K(3.0)
+
+	// Test
+	result := symbolic.VStack(c1, c2, c3, x)
+	if result.Dims()[0] != 4 || result.Dims()[1] != 1 {
+		t.Errorf(
+			"Expected the result to be a vector with 4 elements; received object with shape %v",
+			result.Dims(),
+		)
+	}
+
+	// Verify that the result is a monomial vector
+	if _, ok := result.(symbolic.MonomialVector); !ok {
+		t.Errorf(
+			"Expected the result to be a MonomialVector; received %T",
+			result,
+		)
+	}
+}
+
+/*
+TestExpression_VStack3
+Description:
+
+	Tests the VStack function for 2 vector expressions.
+	Each vector has 11 elements. One is a constant vector and the other is a variable vector.
+	The result should be a monomial vector with 22 elements.
+*/
+func TestExpression_VStack3(t *testing.T) {
+	// Constants
+	kv1 := symbolic.VecDenseToKVector(symbolic.OnesVector(11))
+	vv2 := symbolic.NewVariableVector(11)
+
+	// Test
+	result := symbolic.VStack(kv1, vv2)
+	if result.Dims()[0] != 22 || result.Dims()[1] != 1 {
+		t.Errorf(
+			"Expected the result to be a vector with 22 elements; received object with shape %v",
+			result.Dims(),
+		)
+	}
+
+	// Verify that the result is a monomial vector
+	if _, ok := result.(symbolic.MonomialVector); !ok {
+		t.Errorf(
+			"Expected the result to be a MonomialVector; received %T",
+			result,
+		)
+	}
+}
+
+/*
+TestExpression_VStack4
+Description:
+
+	Tests the VStack function for a matrix and a vector expression.
+	The matrix is a constant matrix and the vector is a variable vector.
+	The matrix is of shape 3x2 and the vector is of length 2.
+	The VStack function should panic because the dimensions do not match.
+*/
+func TestExpression_VStack4(t *testing.T) {
+	// Constants
+	km1 := symbolic.DenseToKMatrix(symbolic.OnesMatrix(3, 2))
+	vv2 := symbolic.NewVariableVector(2)
+
+	// Test
+	defer func() {
+		err := recover().(error)
+		if err == nil {
+			t.Errorf("The VStack function should panic when the dimensions do not match")
+		}
+
+		// Collect the expected error which should be a dimension error and
+		// compare it with the recovered error
+		expectedError := smErrors.DimensionError{
+			Operation: "VStack",
+			Arg1:      km1,
+			Arg2:      vv2,
+		}
+		if err.Error() != expectedError.Error() {
+			t.Errorf(
+				"Expected the error to be %v; received %v",
+				expectedError,
+				err,
+			)
+		}
+
+	}()
+	symbolic.VStack(km1, vv2)
+
+	// The function should panic before this point
+	t.Errorf("The VStack function should panic before this point")
+
+}
+
+/*
+TestExpression_VStack5
+Description:
+
+	Tests the VStack function for a matrix and a matrix expression.
+	The matrix is a constant matrix and the matrix is a variable matrix.
+	The matrix is of shape 3x2 and the matrix is of shape 2x3.
+	The VStack function should panic because the dimensions do not match.
+*/
+func TestExpression_VStack5(t *testing.T) {
+	// Constants
+	km1 := symbolic.DenseToKMatrix(symbolic.OnesMatrix(3, 2))
+	km2 := symbolic.DenseToKMatrix(symbolic.OnesMatrix(2, 3))
+
+	// Test
+	defer func() {
+		err := recover().(error)
+		if err == nil {
+			t.Errorf("The VStack function should panic when the dimensions do not match")
+		}
+
+		// Collect the expected error which should be a dimension error and
+		// compare it with the recovered error
+		expectedError := smErrors.DimensionError{
+			Operation: "VStack",
+			Arg1:      km1,
+			Arg2:      km2,
+		}
+		if err.Error() != expectedError.Error() {
+			t.Errorf(
+				"Expected the error to be %v; received %v",
+				expectedError,
+				err,
+			)
+		}
+
+	}()
+	symbolic.VStack(km1, km2)
+
+	// The function should panic before this point
+	t.Errorf("The VStack function should panic before this point")
+
+}
+
+/*
+TestExpression_VStack6
+Description:
+
+	Tests the VStack function for a matrix and a matrix expression.
+	The matrix is a constant matrix and the matrix is a variable matrix.
+	The matrix is of shape 3x2 and the matrix is of shape 3x2.
+	The VStack function should not panic because the dimensions match.
+	It should produce a monomial matrix with 6 rows and 2 columns.
+*/
+func TestExpression_VStack6(t *testing.T) {
+	// Constants
+	km1 := symbolic.DenseToKMatrix(symbolic.OnesMatrix(3, 2))
+	km2 := symbolic.NewVariableMatrix(3, 2)
+
+	// Test
+	result := symbolic.VStack(km1, km2)
+	if result.Dims()[0] != 6 || result.Dims()[1] != 2 {
+		t.Errorf(
+			"Expected the result to be a 6x2 matrix; received %v",
+			result.Dims(),
+		)
+	}
+
+	// Verify that the result is a monomial matrix
+	if _, ok := result.(symbolic.MonomialMatrix); !ok {
+		t.Errorf(
+			"Expected the result to be a MonomialMatrix; received %T",
+			result,
+		)
+	}
+}