diff --git a/+quantity/Discrete.m b/+quantity/Discrete.m
index fb1663e58164a7a1fe811eb9c8153cb382606b5d..1b832bd1ba73500a5d67728cb874883abf9bbe10 100644
--- a/+quantity/Discrete.m
+++ b/+quantity/Discrete.m
@@ -2151,8 +2151,12 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete ...
 			end
 		end % abs()
 		
-		function d = det(obj)
+		function d = det(obj, setName)
 			% det(X) returns the the determinant of the squre matrix X
+			arguments
+				obj;
+				setName = true;
+			end
 			if isempty(obj)
 				d = quantity.Discrete.empty(1);
 			elseif numel(obj) == 1
@@ -2168,12 +2172,21 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete ...
 						- prod(obj(logical([1, 0, 0; 0, 0, 1; 0, 1, 0]))) ...
 						- prod(obj(logical([0, 1, 0; 1, 0, 0; 0, 0, 1])));
 				else
-					error("only implemented for 2x2 matrices yet");
 					% maybe a more sophisticated implementation would be nice.
+					d = 0;
+					for it = 1 : size(obj, 1)
+						selector = true(size(obj));
+						selector(:, 1) = false;
+						selector(it, :) = false;
+						d = d + (-1)^(it-1)*obj(it, 1) * det(reshape(obj(selector), size(selector)-1), false);
+					end % for it = 1 : size(obj, 1)
 				end
 			else
 				error("det is only defined for quadratic matrices");
 			end
+			if setName
+				d.setName("det(" + obj(1).name + ")");
+			end
 		end % det()
 		
 		function y = real(obj)
diff --git a/+unittests/+quantity/testDiscrete.m b/+unittests/+quantity/testDiscrete.m
index 8a998e75560a30f2c5196ff0b5675c01e0fdbaf7..c304d9ca127095c171c53b7bfec1fe4717385f2a 100644
--- a/+unittests/+quantity/testDiscrete.m
+++ b/+unittests/+quantity/testDiscrete.m
@@ -15,11 +15,17 @@ end
 function testDet(tc)
 z = quantity.Domain("z", linspace(0, 1, 11));
 zDiscrete = quantity.Discrete(z.grid, 'domain', z);
-myMatrix2x2 = [1, 2; 3, 4] * (1 + 0 * zDiscrete);
-myMatrix3x3 = [1, 2, 3; 2, 3, 4; 3, 4, 5] * (1 + 0 * zDiscrete);
+myMatrix2x2 = [1, 2; 3, 4] * zDiscrete;
+myMatrix3x3 = [1, 2, 3; 2, -3, 4; 3, 4, 5] * zDiscrete;
+myMatrix5x5 = [1, 2, 3, 4, 0; ...
+				2, 3, 4, -5, 6; ...
+				3, 4, 5, 6, 7; ...
+				4, -6, 7, 8, 9; ...
+				1, 1, 1, 1, 1] * zDiscrete;
 tc.verifyEqual(zDiscrete.det.on(), zDiscrete.on(), 'AbsTol', 10*eps);
 tc.verifyEqual(myMatrix2x2.det.at(0.5), det(myMatrix2x2.at(0.5)), 'AbsTol', 10*eps);
 tc.verifyEqual(myMatrix3x3.det.at(0.5), det(myMatrix3x3.at(0.5)), 'AbsTol', 10*eps);
+tc.verifyEqual(myMatrix5x5.det.at(0.5), det(myMatrix5x5.at(0.5)), 'AbsTol', 10*eps);
 end % testDet()
 
 function testSum(tc)
diff --git a/+unittests/+quantity/testSymbolic.m b/+unittests/+quantity/testSymbolic.m
index b4e889ae37f27e9d9dcdd131407f60501aab251c..6dbda0d35779a0253beda65e2a6d1882713c6037 100644
--- a/+unittests/+quantity/testSymbolic.m
+++ b/+unittests/+quantity/testSymbolic.m
@@ -7,11 +7,17 @@ end
 function testDet(tc)
 z = quantity.Domain("z", linspace(0, 1, 11));
 zSymbolic = quantity.Symbolic(sym("z"), 'domain', z);
-myMatrix2x2 = [1, 2; 3, 4] * (1 + 0 * zSymbolic);
-myMatrix3x3 = [1, 2, 3; 2, 3, 4; 3, 4, 5] * (1 + 0 * zSymbolic);
+myMatrix2x2 = [1, 2; 3, 4] * zSymbolic;
+myMatrix3x3 = [1, 2, 3; 2, -3, 4; 3, 4, 5] * zSymbolic;
+myMatrix5x5 = [1, 2, 3, 4, 5; ...
+				2, 3, 4, 5, 6; ...
+				3, 4, 5, 6, 7; ...
+				4, -6, 7, 8, 9; ...
+				1, 1, 1, 1, 1] * zSymbolic;
 tc.verifyEqual(zSymbolic.det.on(), zSymbolic.on(), 'AbsTol', 10*eps);
 tc.verifyEqual(myMatrix2x2.det.at(0.5), det(myMatrix2x2.at(0.5)), 'AbsTol', 10*eps);
 tc.verifyEqual(myMatrix3x3.det.at(0.5), det(myMatrix3x3.at(0.5)), 'AbsTol', 10*eps);
+tc.verifyEqual(myMatrix5x5.det.at(0.5), det(myMatrix5x5.at(0.5)), 'AbsTol', 10*eps);
 end % testDet()
 
 function testSum(tc)