diff --git a/+quantity/Discrete.m b/+quantity/Discrete.m
index 650e062a842267a307a26218ebbdd11e7ad18495..c261fd9f5abfad5a51e8fe6f4668fd942b7e6095 100644
--- a/+quantity/Discrete.m
+++ b/+quantity/Discrete.m
@@ -175,11 +175,19 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
 		% --- getter and setters ---
 		%---------------------------
 		function gridName = get.gridName(obj)
-			gridName = {obj.domain.name};
+			if isempty(obj.domain)
+				gridName = {};
+			else
+				gridName = {obj.domain.name};
+			end
 		end
 		
 		function grid = get.grid(obj)
-			grid = {obj.domain.grid};
+			if isempty(obj.domain)
+				grid = {};
+			else
+				grid = {obj.domain.grid};
+			end
 		end
 		
 		function itIs = isConstant(obj)
@@ -190,23 +198,6 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
 		function doNotCopy = get.doNotCopy(obj)
 			doNotCopy = obj.doNotCopyPropertiesName();
 		end
-		function set.gridName(obj, name)
-			if ~iscell(name)
-				name = {name};
-			end
-			obj.gridName = name;
-		end
-		
-		function set.grid(obj, grid)
-			if ~iscell(grid)
-				grid = {grid};
-			end
-			
-			isV = cellfun(@(v) (sum(size(v)>1) == 1) || (numel(v) == 1), grid); % also allow 1x1x1x41 grids
-			assert(all(isV(:)), 'Please use vectors for the grid entries!');
-			
-			[obj.grid] = deal(grid);
-		end
 		function valueDiscrete = get.valueDiscrete(obj)
 			% check if the value discrete for this object
 			% has already been computed.
@@ -430,39 +421,19 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
 			% obj = sortGrid(obj) sorts the grid in alphabetical order.
 			% obj = sort(obj, 'descend') sorts the grid in descending
 			% alphabetical order.
-			
-			if nargin == 2 && strcmp(varargin{1}, 'descend')
-				descend = 1;
-			else
-				descend = 0;
-			end
-			
+						
 			% only sort the grids if there is something to sort
 			if obj(1).nargin > 1
-				gridNames = obj(1).gridName;
-				
-				% this is the default case for ascending alphabetical
-				% order
-				[sortedNames, I] = sort(gridNames);
 				
-				% if descending: flip the order of the entries
-				if descend
-					sortedNames = flip(sortedNames);
-					I = flip(I);
-				end
-				
-				% sort the grid entries
-				[obj.grid] = deal(obj(1).grid(I));
+				[sortedDomain, I] = obj(1).domain.sort(varargin{:});
+				[obj.domain] = deal(sortedDomain);
 				
-				% assign the new grid names
-				[obj.gridName] = deal(sortedNames);
-				
-				% permute the value discrete
 				for k = 1:numel(obj)
 					obj(k).valueDiscrete = permute(obj(k).valueDiscrete, I);
 				end
 			end
 		end% sort()
+		
 		function c = horzcat(a, varargin)
 			%HORZCAT Horizontal concatenation.
 			%   [A B] is the horizontal concatenation of objects A and B
@@ -733,15 +704,15 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
 			myParser.addParameter('AbsTol', 1e-6);
 			myParser.parse(varargin{:});
 			
-			variableGrid = myParser.Results.variableGrid;
+			variableGrid = myParser.Results.variableGrid(:);
 			myGridSize = [numel(variableGrid), ...
 				numel(myParser.Results.initialValueGrid)];
 			
 			% the time (s) vector has to start at 0, to ensure the IC. If
 			% variableGrid does not start with 0, it is separated in
 			% negative and positive parts and later combined again.
-			positiveVariableGrid = [0, variableGrid(variableGrid > 0)];
-			negativeVariableGrid = [0, flip(variableGrid(variableGrid < 0))];
+			positiveVariableGrid = [0; variableGrid(variableGrid > 0)];
+			negativeVariableGrid = [0; flip(variableGrid(variableGrid < 0))];
 			
 			% solve ode for every entry in obj and for every initial value
 			options = odeset('RelTol', myParser.Results.RelTol, 'AbsTol', myParser.Results.AbsTol);
@@ -1006,7 +977,7 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
 		function s = gridSize(obj, myGridName)
 			% GRIDSIZE returns the size of all grid entries.
 			% todo: this should be called gridLength
-			if isempty(obj(1).grid)
+			if isempty(obj(1).domain)
 				s = [];
 			else
 				if nargin == 1
@@ -1449,7 +1420,7 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
 			if isscalar(a)
 				% do the scalar multiplication in a loop
 				for k = 1:numel(b)
-					P(k) = innerMTimes(a, b(k));
+					P(k) = innerMTimes(a, b(k)); %#ok<AGROW>
 				end
 				P = reshape(P, size(b));
 				return
@@ -1470,46 +1441,48 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
 			parameters.name = [a(1).name, ' ', b(1).name];
 			parameters.figureID = a(1).figureID;
 			
-			
 			domainA = a(1).domain;
 			domainB = b(1).domain;
 			
-			% select finest grid from a and b
-			joinedDomain = gridJoin(domainA, domainB);
-			newDomainName = {joinedDomain.name};
-			
-			% gridJoin combines the grid of a and b while using the finer
-			% of both.
-			
-			finalGridOrder = [domainA(idx.A.common), domainA(~idx.A.common), domainB(~idx.B.common)];
+			% The joined domain of quantity A and B is a mixture from both
+			% domains. If the domains have the same name, it must be
+			% checked if they have the same discretization. This is done by
+			% the function "gridJoin". It returns the finer grid.
+			% If the domain names do not coincide, they are just
+			% appended. At first, the domains of qunatity A, then the
+			% domains of quantity B. 
+			joinedDomain = [ gridJoin(domainA(idx.A.common), domainB(idx.B.common)), ...
+				domainA(~idx.A.common), domainB(~idx.B.common) ];
 			
 			%% generate the new grids
 			newDomainA = quantity.Domain.empty(); % for later call of a.on() with fine grid
 			newDomainB = quantity.Domain.empty(); % for later call of b.on() with fine grid
-			for i = 1 : numel(newDomainName)
+			for i = 1 : numel(joinedDomain)
 
-				parameters.domain(i) = finalGridOrder(i);
+				parameters.domain(i) = joinedDomain(i);
 				
 				% if there is a component of the original domain in the
 				% joined domain, it can happen the length of a domain has
 				% changed. Thus, it must be updated for a later evaluation.
-				idxA = domainA.gridIndex(finalGridOrder(i).name);
-				idxB = domainB.gridIndex(finalGridOrder(i).name);
+				idxA = domainA.gridIndex(joinedDomain(i).name);
+				idxB = domainB.gridIndex(joinedDomain(i).name);
 				
 				if idxA
-					newDomainA = [newDomainA, finalGridOrder(i)];
+					newDomainA = [newDomainA, joinedDomain(i)]; %#ok<AGROW>
 				end
 				
 				if idxB
-					newDomainB = [newDomainB, finalGridOrder(i)];
+					newDomainB = [newDomainB, joinedDomain(i)]; %#ok<AGROW>
 				end
 			end
 			parameters = misc.struct2namevaluepair(parameters);
 			
-			%% 
+			% evaluation of the quantities on their "maybe" new domain.
 			valueA = a.on(newDomainA);
 			valueB = b.on(newDomainB);
 			
+			% do the multidimensional tensor multiplication and permute the
+			% values to the right order 
 			C = misc.multArray(valueA, valueB, idx.A.value(end), idx.B.value(1), idx.common);
 			C = permute(C, permuteGrid);
 			P = quantity.Discrete(C, parameters{:});
@@ -1674,11 +1647,13 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
 			% Check if there is any dimension which is zero
 			empty = any(size(obj) == 0);
 			
-			% If the constructor is called without arguments, a
-			% quantity.Discrete is created without an initialization of
-			% the grid. Thus, the quantity is not initialized and empty if
-			% the grid is not a cell.
-			empty = empty || ~iscell(obj(1).grid);
+			% If the constructor is called without any arguments, an object
+			% is still created ( this is required to allow object-arrays),
+			% but the object should be recognized as an empty object. Thus,
+			% we can test if the domain has been set with a domain object.
+			% This happens only in the part of the constructor which is
+			% entered if a valid quantity is initialized.
+			empty = empty || ~isa(obj(1).domain, 'quantity.Domain');
 			
 		end % isempty()
 		
@@ -2148,17 +2123,21 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
 			end
 			value = reshape(v, [obj(1).gridSize(), size(obj)]);
 			
-			if myDomain ~= obj(1).domain
+			if myDomain ~= obj(1).domain 
 				% if a new domain is specified for the evaluation of
-				% the quantity, do an interpolation based on the old
-				% data.
+				% the quantity, ...
+				if obj.isConstant
+					% ... duplicate the constant value on the desired
+					% domain
+					value = repmat(value, [cellfun(@(v) numel(v), {myDomain.grid}), ones(1, length(size(obj)))]);
+				else
+					%... do an interpolation based on the old data.
 				indexGrid = arrayfun(@(s) 1:s, size(obj), 'UniformOutput', false);
 				tempInterpolant = numeric.interpolant(...
 					[obj(1).grid, indexGrid{:}], value);
 				tempGrid = {myDomain.grid};
 				value = tempInterpolant.evaluate(tempGrid{:}, indexGrid{:});
-			elseif obj.isConstant
-				value = repmat(value, [cellfun(@(v) numel(v), {myDomain.grid}), ones(1, length(size(obj)))]);
+				end
 			end
 			
 		end
@@ -2349,10 +2328,6 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
 			% copied.
 		end
 		
-		% 		function s = getHeader(obj)
-		% 			s = 'TEST'
-		% 		end
-		
 		function s = getPropertyGroups(obj)
 			% Function to display the correct values
 			
@@ -2367,7 +2342,7 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
 				s.PropertyList.valueDiscrete = ...
 					[sprintf('%ix', gridSize(obj), size(obj)) sprintf('\b')];
 			end
-		end
+		end % getPropertyGroups
 	end % methods (Access = protected)
 	
 	methods ( Static, Access = protected )
diff --git a/+quantity/Domain.m b/+quantity/Domain.m
index fb8647eed11b52e69942e617b2ff9e9b2af21e63..ef6fd4100a7e77def3bcedd91d68976f5ee69e63 100644
--- a/+quantity/Domain.m
+++ b/+quantity/Domain.m
@@ -1,4 +1,4 @@
-classdef Domain
+classdef Domain < matlab.mixin.CustomDisplay
 	%DOMAIN class to describes a range of values on which a function can be defined.
 	
 	% todo:
@@ -15,8 +15,6 @@ classdef Domain
 		% a speaking name for this domain; Should be unique, so that the
 		% domain can be identified by the name.
 		name char;
-		
-		isConstant = false ;
 	end
 	
 	properties (Dependent)
@@ -24,6 +22,10 @@ classdef Domain
 		lower;		% lower bound of the domain
 		upper;		% upper bound of the domain
 	end
+
+	properties (Dependent, Hidden)
+		ones;
+	end
 	
 	methods
 		function obj = Domain(varargin)
@@ -45,6 +47,10 @@ classdef Domain
 			end
 		end
 		
+		function i = get.ones(obj)
+			i = ones(size(obj.grid));
+		end
+		
 		function n = get.n(obj)
 			n = length(obj.grid);
 		end
@@ -86,12 +92,12 @@ classdef Domain
 			l = ~(obj == obj2);
 		end
 		
-		function [joinedDomain] = gridJoin(obj1, obj2)
+		function [joinedDomain, index] = gridJoin(obj1, obj2)
 			%% gridJoin combines the grid and gridName of two objects (obj1,
 			% obj2), such that every gridName only occurs once and that the
 			% finer grid of both is used.
 			
-			[joinedGrid] = union({obj1.name}, {obj2.name}, 'sorted');
+			[joinedGrid, index] = unique({obj1.name, obj2.name}, 'stable');
 			
 			joinedDomain = quantity.Domain.empty;
 			
@@ -115,10 +121,9 @@ classdef Domain
 					tempDomain1 = obj1(index1);
 					tempDomain2 = obj2(index2);
 					
-					if ~tempDomain1.isConstant && ~tempDomain2.isConstant
-						assert(tempDomain1.lower == tempDomain2.lower, 'Grids must have same domain for gridJoin')
-						assert(tempDomain1.upper == tempDomain2.upper, 'Grids must have same domain for gridJoin')
-					end
+					assert(tempDomain1.lower == tempDomain2.lower, 'Grids must have same domain for gridJoin')
+					assert(tempDomain1.upper == tempDomain2.upper, 'Grids must have same domain for gridJoin')
+
 					if tempDomain1.gridLength > tempDomain2.gridLength
 						joinedDomain(i) = tempDomain1;
 					else
@@ -192,8 +197,83 @@ classdef Domain
 				newDomain = obj(idx);
 			end
 		end
+		
+		
+		function [newObj, I] = sort(obj, varargin)
+			%SORT sorts the grid of the object in a desired order
+			% obj = sortGrid(obj) sorts the grid in alphabetical order.
+			% obj = sort(obj, 'descend') sorts the grid in descending
+			% alphabetical order.
+			
+			if nargin == 2 && strcmp(varargin{1}, 'descend')
+				descend = 1;
+			elseif nargin == 1 || strcmp(varargin{1}, 'ascend')
+				descend = 0;
+			else
+				error(['Unknown sorting order' char( varargin{1} )])
+			end
+			
+			% only sort the grids if there is something to sort
+			if obj.ndims > 1
+				gridNames = {obj.name};
+				
+				% this is the default case for ascending alphabetical
+				% order
+				[sortedNames, I] = sort(gridNames);
+				
+				% if descending: flip the order of the entries
+				if descend
+					sortedNames = flip(sortedNames);
+					I = flip(I);
+				end
+				
+				% sort the grid entries
+				myGrids = {obj.grid};
+				myGrids = myGrids(I);
+
+				for k = 1:length(myGrids)
+					newObj(k) = quantity.Domain('grid', myGrids{k}, ...
+						'name', sortedNames{k}); %#ok<AGROW>
+				end
+			else
+				newObj = obj;
+			end
+		end% sort()			
 	end
 	
+	methods (Access = protected)
+		
+		function s = getPropertyGroups(obj)
+			% Function to display the correct values
+			
+			if isempty(obj)
+				s = getPropertyGroups@matlab.mixin.CustomDisplay(obj);
+				return;
+			else
+				s = getPropertyGroups@matlab.mixin.CustomDisplay(obj(1));
+			end
+			
+			if numel(obj) ~= 1
+				
+				grids = {obj.grid};
+				
+				sizes = cellfun(@(g)size(g), grids, 'UniformOutput', false);							
+								
+				s.PropertyList.grid = ...
+					[sprintf('%ix%i,  ', cell2mat(sizes)) sprintf('\b\b\b')];
+				s.PropertyList.name = ...
+					[sprintf('%s,  ', obj.name) sprintf('\b\b\b')];
+				s.PropertyList.lower = ... 
+					[sprintf('%d,  ', obj.lower) sprintf('\b\b\b')];
+				s.PropertyList.upper = ... 
+					[sprintf('%.3f,  ', obj.upper) sprintf('\b\b\b')];
+				s.PropertyList.n = ...
+					[sprintf('%i,  ', obj.n) sprintf('\b\b\b')];
+			end
+		end % getPropertyGroups
+		
+	end % Access = protected
+	
 	methods (Static)
 		function g = defaultGrid(gridSize, name)
 			
diff --git a/+quantity/Symbolic.m b/+quantity/Symbolic.m
index f65d413b5fa1393f0a1e10f99719deb5f60a8e09..ba4b12b07f0ec1a9adbda9e893b837658e1b26d5 100644
--- a/+quantity/Symbolic.m
+++ b/+quantity/Symbolic.m
@@ -576,28 +576,28 @@ classdef Symbolic < quantity.Function
 			mObj = copy(obj);
 		end % uplus()
 		
-		function parameters = combineGridGridNameVariable(A, B, parameters)
-			assert(isa(A, 'quantity.Symbolic') && isa(B, 'quantity.Symbolic'), ...
-				'A and B must be quantity.Symbolic');
-			
-			if nargin < 3
-				parameters = struct(A(1));
-			end
-			% combine grids of two quantities. This is often used in mathematical
-			% operations to obtain the parameters of the resulting quantity.
-			idx = computePermutationVectors(A, B);
-			parameters.grid = ...
-				[A(1).grid(idx.A.common), ...
-				A(1).grid(~idx.A.common), ...
-				B(1).grid(~idx.B.common)];
-			parameters.gridName = [A(1).gridName(idx.A.common), ...
-				A(1).gridName(~idx.A.common), ...
-				B(1).gridName(~idx.B.common)];
-			parameters.variable = [A(1).variable(idx.A.common), ...
-				A(1).variable(~idx.A.common), ...
-				B(1).variable(~idx.B.common)];
-		end % combineGridGridNameVariable()
-		
+% 		function parameters = combineGridGridNameVariable(A, B, parameters)
+% 			assert(isa(A, 'quantity.Symbolic') && isa(B, 'quantity.Symbolic'), ...
+% 				'A and B must be quantity.Symbolic');
+% 			
+% 			if nargin < 3
+% 				parameters = struct(A(1));
+% 			end
+% 			% combine grids of two quantities. This is often used in mathematical
+% 			% operations to obtain the parameters of the resulting quantity.
+% 			idx = computePermutationVectors(A, B);
+% 			parameters.grid = ...
+% 				[A(1).grid(idx.A.common), ...
+% 				A(1).grid(~idx.A.common), ...
+% 				B(1).grid(~idx.B.common)];
+% 			parameters.gridName = [A(1).gridName(idx.A.common), ...
+% 				A(1).gridName(~idx.A.common), ...
+% 				B(1).gridName(~idx.B.common)];
+% 			parameters.variable = [A(1).variable(idx.A.common), ...
+% 				A(1).variable(~idx.A.common), ...
+% 				B(1).variable(~idx.B.common)];
+% 		end % combineGridGridNameVariable()
+% 		
 		function C = mtimes(A, B)
 			% if one input is ordinary matrix, this is very simple.
 			if isempty(B) ||  isempty(A)
@@ -624,7 +624,7 @@ classdef Symbolic < quantity.Function
 				return;
 			end
 			
-			parameters = combineGridGridNameVariable(A, B);
+			parameters.domain = gridJoin(A.domain, B.domain);
 			parameters.name = [A(1).name, ' ', B(1).name];
 			parameters = misc.struct2namevaluepair(parameters);
 			
@@ -644,7 +644,7 @@ classdef Symbolic < quantity.Function
 				parameters = struct(a(1));
 			end
 			
-			parameters = combineGridGridNameVariable(a, b, parameters);
+			parameters.domain = gridJoin(a.domain, b.domain);
 			parameters.name = [a(1).name, '+' , b(1).name];
 			parameters = misc.struct2namevaluepair(parameters);
 			
diff --git a/+unittests/+quantity/testDiscrete.m b/+unittests/+quantity/testDiscrete.m
index e18b0fc973d26609401aaa282130e7fbd8fa5031..66da1d4d827123cf04bc1e47df3d718bd2b05bc0 100644
--- a/+unittests/+quantity/testDiscrete.m
+++ b/+unittests/+quantity/testDiscrete.m
@@ -24,12 +24,37 @@ tc.verifyEqual(blubQuadraticNorm.on(), 2*ones(11,1));
 tc.verifyEqual(blubQuadraticWeightedNorm.on(), 4*ones(11,1));
 end % testQuadraticNorm()
 
+function testSort(tc)
+
+t = linspace(0, pi, 7)';
+domA = quantity.Domain('grid', t, 'name', 'a');
+domB = quantity.Domain('grid', t, 'name', 'b');
+q = quantity.Discrete(sin(t) .* cos(t'), ...
+	'domain', [domA, domB]);
+
+q1V = q.valueDiscrete();
+
+q.sort('descend');
+
+q2V = q.valueDiscrete();
+
+tc.verifyEqual(q1V , q2V')
+
+
+end % testSort
+
 function testBlkdiag(tc)
 % init some data
 syms z zeta
-A = quantity.Symbolic(...
-	[1+z*zeta, -zeta; -z, z^2], 'grid', {linspace(0, 1, 21), linspace(0, 1, 41)},...
-	'variable', {z, zeta}, 'name', 'q');
+
+z = quantity.Domain('grid', linspace(0, 1, 21), 'name', 'z');
+zeta = quantity.Domain('grid', linspace(0, 1, 41), 'name', 'zeta');
+
+A = quantity.Discrete(cat(4, ...
+	cat(3, 1 + z.grid .* zeta.grid', - z.ones .* zeta.grid'), ...
+	cat(3, -z.grid .* zeta.ones', z.grid .^2 .* zeta.ones') ), ...
+	'domain', [z, zeta], 'name', 'q');
+
 B = 2*A(:, 1);
 C = 3*A(1);
 
@@ -82,12 +107,14 @@ tc.verifyEqual(qDiscrete2.imag.on(), -qDiscrete2Ctransp.imag.on());
 end % testCtranspose
 
 function testTranspose(tc)
+
 syms z zeta
 qSymbolic = quantity.Symbolic(...
 	[1+z*zeta, -zeta; -z, z^2], 'grid', {linspace(0, 1, 21), linspace(0, 1, 41)},...
 	'variable', {z, zeta}, 'name', 'q');
 qDiscrete = quantity.Discrete(qSymbolic);
 qDiscreteTransp = qDiscrete.';
+
 tc.verifyEqual(qDiscrete(1,1).on(), qDiscreteTransp(1,1).on());
 tc.verifyEqual(qDiscrete(2,2).on(), qDiscreteTransp(2,2).on());
 tc.verifyEqual(qDiscrete(1,2).on(), qDiscreteTransp(2,1).on());
@@ -115,19 +142,22 @@ tc.verifyEqual(fReference2.on(), fFlipped3.on(), 'AbsTol', 10*eps);
 end % testFlipGrid();
 
 function testScalarPlusMinusQuantity(testCase)
-syms z
-myGrid = linspace(0, 1, 7);
-f = quantity.Discrete(quantity.Symbolic([1; 2] + zeros(2, 1)*z, ...
-	'variable', {z}, 'grid', {myGrid}));
+
+z = quantity.Domain('grid', linspace(0, 1, 7), 'name', 'z');
+f = quantity.Discrete( ( [2; 1] + zeros(2, 1) .* z.grid' )', ...
+	'domain', z);
+
 testCase.verifyError(@() 1-f-1, '');
 testCase.verifyError(@() 1+f+1, '');
+
 end % testScalarPlusMinusQuantity
 
 function testNumericVectorPlusMinusQuantity(testCase)
-syms z
-myGrid = linspace(0, 1, 7);
-f = quantity.Discrete(quantity.Symbolic([1+z; 2+sin(z)] + zeros(2, 1)*z, ...
-	'variable', {z}, 'grid', {myGrid}));
+
+z = quantity.Domain('grid', linspace(0, 1, 7), 'name', 'z');
+f = quantity.Discrete( [1+z.grid, 2+sin(z.grid)] , ...
+	'domain', z);
+
 a = ones(size(f));
 testCase.verifyEqual(on(a-f), 1-f.on());
 testCase.verifyEqual(on(f-a), f.on()-1);
@@ -136,19 +166,18 @@ testCase.verifyEqual(on(f+a), f.on()+1);
 end % testNumericVectorPlusMinusQuantity
 
 function testUnitaryPluasAndMinus(testCase)
-syms z zeta
-qSymbolic = quantity.Symbolic(...
-	[1+z*zeta, -zeta; -z, z^2], 'grid', {linspace(0, 1, 21), linspace(0, 1, 41)},...
-	'variable', {z, zeta}, 'name', 'q');
-qDiscrete = quantity.Discrete(qSymbolic);
-qDoubleArray = qSymbolic.on();
+z = quantity.Domain('grid', linspace(0, 1, 7), 'name', 'z');
+zeta = quantity.Domain('grid', linspace(0, 1, 7), 'name', 'zeta');
+
+q = quantity.Discrete(cat(4, ...
+	cat(3, 1 + z.grid .* zeta.grid', - z.grid .* zeta.ones'), ...
+	cat(3, - z.ones .* zeta.grid', z.grid.^2 .* zeta.ones') ), ...
+	'domain', [z, zeta]) ;
 
-testCase.verifyEqual(on(-qSymbolic), -qDoubleArray);
-testCase.verifyEqual(on(-qDiscrete), -qDoubleArray);
-testCase.verifyEqual(on(+qSymbolic), +qDoubleArray);
-testCase.verifyEqual(on(+qDiscrete), +qDoubleArray);
-testCase.verifyEqual(on(+qDiscrete), on(+qSymbolic));
-testCase.verifyEqual(on(-qDiscrete), on(-qSymbolic));
+qDoubleArray = q.on();
+	
+testCase.verifyEqual(on(-q), -qDoubleArray);
+testCase.verifyEqual(on(+q), +qDoubleArray);
 
 end
 
@@ -198,27 +227,33 @@ end
 
 function testExp(testCase)
 % 1 spatial variable
-syms z zeta
-s1d = quantity.Discrete(quantity.Symbolic(...
-	1+z*z, 'grid', {linspace(0, 1, 21)}, 'variable', {z}, 'name', 's1d'));
+z = quantity.Domain('grid', linspace(0, 1, 3), 'name', 'z');
+s1d = quantity.Discrete(1 + z.grid .* z.grid.', 'domain', z);
 testCase.verifyEqual(s1d.exp.on(), exp(s1d.on()));
 
 % diagonal matrix
-s2dDiag = quantity.Discrete(quantity.Symbolic(...
-	[1+z*zeta, 0; 0, z^2], 'grid', {linspace(0, 1, 21), linspace(0, 1, 41)},...
-	'variable', {z, zeta}, 'name', 's2dDiag'));
+zeta = quantity.Domain('grid', linspace(0, 1, 4), 'name', 'zeta');
+s2dDiag = quantity.Discrete(cat(4, ...
+	cat(3, 1 + z.grid .* zeta.grid', zeros(z.n, zeta.n)), ....
+	cat(3, zeros(z.n, zeta.n), z.grid.^2 .* zeta.ones')), 'domain', [z, zeta]);
+
 testCase.verifyEqual(s2dDiag.exp.on(), exp(s2dDiag.on()));
 end
 
 function testExpm(testCase)
-syms z zeta
-mat2d = quantity.Discrete(quantity.Symbolic(...
-	ones(2, 2) + [1+z*zeta, 3*zeta; 2+5*z+zeta, z^2], 'grid', {linspace(0, 1, 21), linspace(0, 1, 41)},...
-	'variable', {z, zeta}, 'name', 's2d'));
+
+z = quantity.Domain('grid', linspace(0, 1, 3), 'name', 'z');
+zeta = quantity.Domain('grid', linspace(0, 1, 4), 'name', 'zeta');
+
+mat2d = quantity.Discrete(...
+	{1 + z.grid .* zeta.grid', 3 * z.ones .* zeta.grid'; ...
+	 2 + 5 * z.grid + zeta.grid', z.grid .^2 .* zeta.ones'}, ...
+	 'domain', [z, zeta]);
+
 mat2dMat = mat2d.on();
 mat2dExpm = 0 * mat2d.on();
-for zIdx = 1 : 21
-	for zetaIdx = 1 : 41
+for zIdx = 1 : z.n
+	for zetaIdx = 1 : zeta.n
 		mat2dExpm(zIdx, zetaIdx, :, :) = expm(reshape(mat2dMat(zIdx, zetaIdx, :, :), [2, 2]));
 	end
 end
@@ -302,6 +337,7 @@ testCase.verifyEqual(myMatrixDiscrete(1,1).on(), myVectorDiscrete(1,1).on());
 testCase.verifyEqual(myMatrixDiscrete(2,2).on(), myVectorDiscrete(2,1).on());
 testCase.verifyEqual(numel(myVectorDiscrete), size(myMatrixDiscrete, 1));
 end
+
 function testVec2Diag(testCase)
 % quantity.Discrete
 n = 7;
@@ -343,39 +379,39 @@ end
 
 function testSolveDVariableEqualQuantityConstant(testCase)
 %% simple constant case
-quan = quantity.Discrete(ones(51, 1), 'grid', linspace(0, 1, 51), ...
+quan = quantity.Discrete(ones(3, 1), 'grid', linspace(0, 1, 3), ...
 	'gridName', 'z', 'size', [1, 1], 'name', 'a');
 solution = quan.solveDVariableEqualQuantity();
-[referenceResult1, referenceResult2] = ndgrid(linspace(0, 1, 51), linspace(0, 1, 51));
+[referenceResult1, referenceResult2] = ndgrid(linspace(0, 1, 3), linspace(0, 1, 3));
 testCase.verifyEqual(solution.on(), referenceResult1 + referenceResult2, 'AbsTol', 10*eps);
 end
 
 function testSolveDVariableEqualQuantityNegative(testCase)
 syms z zeta
 assume(z>0 & z<1); assume(zeta>0 & zeta<1);
-myParameterGrid = linspace(0, 1, 51);
+myParameterGrid = linspace(0, 0.1, 5);
 Lambda = quantity.Symbolic(-0.1-z^2, ...%, -1.2+z^2]),...1+z*sin(z)
 	'variable', z, ...
 	'grid', myParameterGrid, 'gridName', 'z', 'name', '\Lambda');
 
 %%
-myGrid = linspace(-2, 2, 101);
-sEval = -0.9;
+myGrid = linspace(-.2, .2, 11);
+sEval = -0.09;
 solveThisOdeDiscrete = solveDVariableEqualQuantity(...
 	diag2vec(quantity.Discrete(Lambda)), 'variableGrid', myGrid);
 solveThisOdeSymbolic = solveDVariableEqualQuantity(...
 	diag2vec(Lambda), 'variableGrid', myGrid);
 
-testCase.verifyEqual(solveThisOdeSymbolic.on({sEval, 0.5}), solveThisOdeDiscrete.on({sEval, 0.5}), 'RelTol', 5e-4);
+testCase.verifyEqual(solveThisOdeSymbolic.on({sEval, 0.05}), solveThisOdeDiscrete.on({sEval, 0.05}), 'RelTol', 5e-3);
 end
 
 function testSolveDVariableEqualQuantityComparedToSym(testCase)
 %% compare with symbolic implementation
 syms z
 assume(z>0 & z<1);
-quanBSym = quantity.Symbolic(1+z, 'grid', {linspace(0, 1, 21)}, ...
+quanBSym = quantity.Symbolic(1+z, 'grid', {linspace(0, 1, 5)}, ...
 	'gridName', 'z', 'name', 'bSym', 'variable', {z});
-quanBDiscrete = quantity.Discrete(quanBSym.on(), 'grid', {linspace(0, 1, 21)}, ...
+quanBDiscrete = quantity.Discrete(quanBSym.on(), 'grid', {linspace(0, 1, 5)}, ...
 	'gridName', 'z', 'name', 'bDiscrete', 'size', size(quanBSym));
 solutionBSym = quanBSym.solveDVariableEqualQuantity();
 solutionBDiscrete = quanBDiscrete.solveDVariableEqualQuantity();
@@ -387,10 +423,11 @@ function testSolveDVariableEqualQuantityAbsolut(testCase)
 %% compare with symbolic implementation
 syms z
 assume(z>0 & z<1);
-quanBSym = quantity.Symbolic(1+z, 'grid', {linspace(0, 1, 51)}, ...
+quanBSym = quantity.Symbolic(1+z, 'grid', {linspace(0, 1, 11)}, ...
 	'gridName', 'z', 'name', 'bSym', 'variable', {z});
-quanBDiscrete = quantity.Discrete(quanBSym.on(), 'grid', {linspace(0, 1, 51)}, ...
+quanBDiscrete = quantity.Discrete(quanBSym.on(), 'grid', {linspace(0, 1, 11)}, ...
 	'gridName', 'z', 'name', 'bDiscrete', 'size', size(quanBSym));
+
 solutionBDiscrete = quanBDiscrete.solveDVariableEqualQuantity();
 myGrid = solutionBDiscrete.grid{1};
 solutionBDiscreteDiff = solutionBDiscrete.diff(1, 's');
@@ -399,9 +436,10 @@ for icIdx = 1 : length(myGrid)
 	quanOfSolutionOfS(:, icIdx, :) = quanBSym.on(solutionBDiscrete.on({myGrid, myGrid(icIdx)}));
 end
 %%
-testCase.verifyEqual(solutionBDiscreteDiff.on({myGrid(2:end-1), myGrid}), quanOfSolutionOfS(2:end-1, :), 'AbsTol', 1e-3);
+testCase.verifyEqual(solutionBDiscreteDiff.on({myGrid(2:end-1), myGrid}), quanOfSolutionOfS(2:end-1, :), 'AbsTol', 1e-2);
 
 end
+
 function testMtimesDifferentGridLength(testCase)
 %%
 a = quantity.Discrete(ones(11, 1), 'grid', linspace(0, 1, 11), 'gridName', 'z', ...
@@ -414,7 +452,6 @@ syms z
 c = quantity.Symbolic(1, 'grid', linspace(0, 1, 21), 'variable', z, 'name', 'c');
 ac = a*c;
 
-%%
 %%
 testCase.verifyEqual(ab.on(), ones(11, 1));
 testCase.verifyEqual(ac.on(), ones(21, 1));
@@ -523,17 +560,24 @@ testCase.verifyEqual([a(1).on({zeta(2), eta(3), z(4)}); a(2).on({zeta(2), eta(3)
 end
 
 function testMtimesZZeta2x2(testCase)
-gridVecA = linspace(0, 1, 101);
-[zGridA , ~] = ndgrid(gridVecA, gridVecA);
 
-a = quantity.Discrete(ones([size(zGridA), 2, 2]), 'size', [2, 2], ...
-	'grid', {gridVecA, gridVecA}, 'gridName', {'z', 'zeta'}, 'name', 'A');
+
+
+%%
+z11_ = linspace(0, 1, 11);
+z11 = quantity.Domain('grid', z11_, 'name', 'z');
+zeta11 = quantity.Domain('grid', z11_, 'name', 'zeta');
+
+a = quantity.Discrete(ones([z11.n, z11.n, 2, 2]), 'size', [2, 2], ...
+	'domain', [z11, zeta11], 'name', 'A');
 
 syms zeta
-assume(zeta>=0 & zeta>=1)
-gridVecB = linspace(0, 1, 100);
+z10_ = linspace(0, 1, 10);
+z10 = quantity.Domain('grid', z10_, 'name', 'zeta');
+
 b = quantity.Symbolic((eye(2, 2)), 'variable', zeta, ...
-	'grid', gridVecB);
+	'grid', z10_);
+
 c = a*b;
 
 %
@@ -543,10 +587,11 @@ end
 function testMTimesPointWise(tc)
 
 syms z zeta
-Z = linspace(0, 1, 501)';
-ZETA = Z;
-P = quantity.Symbolic([z, z^2; z^3, z^4] * zeta, 'grid', {Z, ZETA});
-B = quantity.Symbolic([sin(zeta); cos(zeta)], 'grid', ZETA);
+Z = linspace(0, pi, 51)';
+ZETA = linspace(0, pi / 2, 71)';
+
+P = quantity.Symbolic( sin(zeta) * z, 'grid', {Z, ZETA}, 'gridName', {'z', 'zeta'})
+B = quantity.Symbolic( sin(zeta), 'grid', ZETA, 'gridName', 'zeta');
 
 PB = P*B;
 
@@ -554,6 +599,7 @@ p = quantity.Discrete(P);
 b = quantity.Discrete(B);
 
 pb = p*b;
+
 tc.verifyEqual(MAX(abs(PB-pb)), 0, 'AbsTol', 10*eps);
 end
 
diff --git a/+unittests/+quantity/testDomain.m b/+unittests/+quantity/testDomain.m
index da316cad0b01cbea9e3daebea28a35ec9cf86b63..ea22a6db4a81c48e3d6104522dd3e71535e4376d 100644
--- a/+unittests/+quantity/testDomain.m
+++ b/+unittests/+quantity/testDomain.m
@@ -15,6 +15,26 @@ testCase.verifyEqual( ndims(D), 2);
 testCase.verifyEqual( ndims(d), 1);
 end
 
+function testSort(testCase)
+
+z = linspace(0, pi, 3);
+a = quantity.Domain('grid', z, 'name', 'a');
+b = quantity.Domain('grid', z, 'name', 'b');
+c = quantity.Domain('grid', z, 'name', 'c');
+
+d = [a b c];
+
+[d_, I_] = d.sort('descend');
+[d__, I__] = d.sort('ascend');
+d___ = d.sort();
+
+testCase.verifyEqual({d_.name}, {'c', 'b', 'a'});
+testCase.verifyEqual({d__.name}, {d.name});
+testCase.verifyEqual({d__.name}, {d___.name});
+testCase.verifyEqual( I_, flip(I__) );
+
+
+end
 
 function testDomainGridLength(testCase)
 
@@ -75,9 +95,10 @@ c = Pt;
 % c = quantity.Discrete(ones(numel(t), 2, 2), ...
 % 	'size', [2 2], 'grid', {t}, 'gridName', {'p'});
 
-joinedDomainAB = gridJoin(a, b);
-joinedDomainCC = gridJoin(c, c);
-joinedDomainBC = gridJoin(b, c);
+joinedDomainAB = sort( gridJoin(a, b) );
+joinedDomainCC = sort( gridJoin(c, c) );
+joinedDomainBC = sort( gridJoin(b, c) );
+joinedDomain = sort( gridJoin(Z, [Z, T]) );
 
 testCase.verifyEqual({joinedDomainAB.name}, {'p', 's', 't', 'z'});
 testCase.verifyEqual({joinedDomainAB.grid}, {s(:), s(:), t(:), z(:)});
@@ -85,6 +106,7 @@ testCase.verifyEqual({joinedDomainCC.name}, {'p'});
 testCase.verifyEqual({joinedDomainCC.grid}, {t(:)});
 testCase.verifyEqual({joinedDomainBC.name}, {'p', 's'});
 testCase.verifyEqual({joinedDomainBC.grid}, {s(:), t(:)});
+testCase.verifyEqual({joinedDomain.grid}, {t(:), z(:)});
 end
 
 function testGridIndex(testCase)
diff --git a/+unittests/+quantity/testSymbolic.m b/+unittests/+quantity/testSymbolic.m
index 5634c445b328a6105b88f0d24e20ed1299677aa2..dc68188de4f00cc74656588918ce5f5225165f89 100644
--- a/+unittests/+quantity/testSymbolic.m
+++ b/+unittests/+quantity/testSymbolic.m
@@ -611,20 +611,26 @@ verifyEqual(testCase, q.variable, z);
 
 end
 
-function testUMinus(testCase)
+function testUMinusAndUPlus(testCase)
 %%
 syms z;
 f1 = (sin(z .* pi) );
 f2 =  z;
 f3 = (cos(z .* pi) );
 zz = linspace(0, 1, 42).';
+
 f = quantity.Symbolic([f1, f2; 0, f3], ...
 	'grid', {zz});
-
-mf = -f;
+fDisc = quantity.Discrete(f);
+values = f.on();
 
 %%
-verifyEqual(testCase, [f.valueDiscrete], -[mf.valueDiscrete]);
+testCase.verifyEqual(on( -f ), -values);
+testCase.verifyEqual(on( + ( - f ) ), values );
+testCase.verifyEqual(on( + fDisc), on(+f));
+testCase.verifyEqual(on( - fDisc), on(-f));
+
+
 end
 
 function testConstantValues(testCase)