diff --git a/+quantity/Discrete.m b/+quantity/Discrete.m
index b1ae25f139029fc2be293953f9be47e6b70474b1..96490e08ed98ccc98c353cf86b169239defc61ce 100644
--- a/+quantity/Discrete.m
+++ b/+quantity/Discrete.m
@@ -1609,24 +1609,10 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete ...
 		end % sqrtm()
 		
 		function P = power(obj, p)
-			% a.^p elementwise power implemented by multiplication
-			% #todo unittest required
-			P = quantity.Discrete.zeros(size(obj), obj(1).domain);
-			
-			if round(p) == p
-				% power of integer order:
-				for k = 1:numel(obj)
-					P(k) = obj(k)^p;
-				end
-			else
-				% power of not integer order
-				for k = 1:numel(obj)
-					P(k) = obj(k).copy();
-					P(k).valueDiscrete = obj(k).valueDiscrete.^p;
-					P(k).name = obj(k).name + ".^" + num2str(p);
-				end
-			end
-		end % power 
+			% a.^p elementwise power
+			P = quantity.Discrete(obj.on().^(p), obj(1).domain, ...
+				"name", obj(1).name + ".^{" + num2str(p) + "}");
+		end % power()
 		
 		function P = mpower(a, p)
 			% Matrix power a^p is matrix or scalar a to the power p.
diff --git a/+unittests/+quantity/testDiscrete.m b/+unittests/+quantity/testDiscrete.m
index 1a1014ff8a87e05d309dd7b096f753c868d4a592..1cf32c1c76f9eb12d0e0d382cd2964bb17eafc5a 100644
--- a/+unittests/+quantity/testDiscrete.m
+++ b/+unittests/+quantity/testDiscrete.m
@@ -4,6 +4,29 @@ function [tests] = testDiscrete()
 tests = functiontests(localfunctions);
 end
 
+function testPower(tc)
+%% scalar case
+z = quantity.Domain("z", linspace(0, 1, 11));
+A = quantity.Discrete(sin(z.grid * pi), z);
+aa = sin(z.grid * pi) .* sin(z.grid * pi);
+AA = A.^2;
+
+tc.verifyEqual( aa, AA.on());
+tc.verifyEqual( on(A.^0), ones(z.n, 1));
+tc.verifyEqual( on(A.^3), on(AA * A), 'AbsTol', 1e-15);
+
+%% matrix case
+zeta = quantity.Domain("zeta", linspace(0, 1, 21));
+B = [0*z.Discrete + zeta.Discrete, A + 0*zeta.Discrete; ...
+	z.Discrete+zeta.Discrete, A * subs(A, "z", "zeta")];
+BBB = B.^3;
+BBB_alt = B .* B .* B;
+tc.verifyEqual(BBB.on(), BBB_alt.on(), 'AbsTol', 1e-14);
+tc.verifyEqual(BBB.at({0.5, 0.6}), (B.at({0.5, 0.6})).^3, 'AbsTol', 1e-15);
+tc.verifyEqual(on(B.^0), ones([B(1).domain(1).n, B(1).domain(2).n, 2, 2]), ...
+	'AbsTol', 1e-15);
+end % testPower()
+
 function setupOnce(testCase)
 syms z zeta t sy
 testCase.TestData.sym.z = z;
diff --git a/+unittests/+quantity/testSymbolic.m b/+unittests/+quantity/testSymbolic.m
index 52066431560ecaf8fef302f2aee53ead3b946f7a..25f1b37226771476fac2ff023dfe6ae4754010f0 100644
--- a/+unittests/+quantity/testSymbolic.m
+++ b/+unittests/+quantity/testSymbolic.m
@@ -13,12 +13,35 @@ B = [0*z.Symbolic + zeta.Symbolic, A + 0*zeta.Symbolic; ...
 B2 = [zeta.Symbolic + 0*z.Symbolic, A + 0*zeta.Symbolic; ...
 	z.Symbolic+zeta.Symbolic, A * subs(A, "z", "zeta")];
 tc.verifyEqual(B.on([z, zeta]), B2.on([z, zeta]));
-% TODO: fix this error!
+% TODO: fix this error! See issue #41
 % this problem is maybe caused as input arguments of B2(1).valueContinuous is somehow ordered
 % differently compared to B2(2:4).valueContinuous or all B(:).valueContinuous, 
 % i.e. B2(1).valueContinuous is (zeta, z) instead of (z, zeta).
 end % testValueContinuousBug()
 
+function testPower(tc)
+%% scalar case
+z = quantity.Domain("z", linspace(0, 1, 11));
+A = quantity.Symbolic(sin(sym("z") * pi), z);
+aa = sin(z.grid * pi) .* sin(z.grid * pi);
+AA = A.^2;
+
+tc.verifyEqual( aa, AA.on());
+tc.verifyEqual( on(A.^0), ones(z.n, 1));
+tc.verifyEqual( on(A.^3), on(AA * A), 'AbsTol', 1e-15);
+
+%% matrix case
+zeta = quantity.Domain("zeta", linspace(0, 1, 21));
+B = [0*z.Symbolic + zeta.Symbolic, A + 0*zeta.Symbolic; ...
+	z.Symbolic+zeta.Symbolic, A * subs(A, "z", "zeta")];
+BBB = B.^3;
+BBB_alt = B .* B .* B;
+tc.verifyEqual(BBB.on(), BBB_alt.on(), 'AbsTol', 1e-14);
+tc.verifyEqual(BBB.at({0.5, 0.6}), (B.at({0.5, 0.6})).^3, 'AbsTol', 1e-15);
+tc.verifyEqual(on(B.^0), ones([B(1).domain(1).n, B(1).domain(2).n, 2, 2]), ...
+	'AbsTol', 1e-15);
+end % testPower()
+
 function testMPower(tc)
 %% scalar case
 z = quantity.Domain("z", linspace(0, 1, 11));