diff --git a/+quantity/Discrete.m b/+quantity/Discrete.m
index 65283338a6f431340cca2cc692ac74a6ace5fcce..e25c98cb6adfde9b940b2acdc12b10859d245fbe 100644
--- a/+quantity/Discrete.m
+++ b/+quantity/Discrete.m
@@ -2158,7 +2158,7 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
 			
 			value = reshape(v, [obj(1).gridSize(), size(obj)]);
 				
-			if nargin == 2 && myDomain ~= obj(1).domain 
+			if nargin == 2 && ~isequal( myDomain, obj(1).domain )
 				% if a new domain is specified for the evaluation of
 				% the quantity, ...
 				if obj.isConstant
@@ -2210,15 +2210,17 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
 		
 		function P = zeros(valueSize, grid, varargin)
 			%ZEROS initializes an zero quantity.Discrete object
-			%   P = zeros(gridSize, valueSize, grid) returns a
-			%   quantity.Discrete object that has only zero entries on a
-			%   grid with defined by the cell GRID and the value size
-			%   defined by the size-vector VALUESIZE
+			%   P = zeros(VALUESIZE, GRID) returns a quantity.Discrete
+			%   object that has only zero entries on a grid which is
+			%   defined by the cell GRID and the value size defined by the
+			%   size-vector VALUESIZE
 			%       P = quantity.Discrete([n,m], {linspace(0,1)',
 			%                                       linspace(0,10)});
 			%       creates an (n times m) zero quantity.Discrete on the
 			%       grid (0,1) x (0,10)
 			%
+			%	P = zeros(VALUESIZE, DOMAIN) creates a matrix of size
+			%	VALUESIZE on the DOMAIN with zero entries.
 			if ~iscell(grid)
 				grid = {grid};
 			end
diff --git a/+quantity/Domain.m b/+quantity/Domain.m
index e8832291056b4e9cf4a4380c5e97cc36c85218e8..f17d5a02c548183420f15dc132491c3f7d3850c3 100644
--- a/+quantity/Domain.m
+++ b/+quantity/Domain.m
@@ -1,4 +1,4 @@
-classdef Domain < matlab.mixin.CustomDisplay
+classdef Domain < handle & matlab.mixin.CustomDisplay 
 	%DOMAIN class to describes a range of values on which a function can be defined.
 	
 	% todo:
@@ -42,11 +42,10 @@ classdef Domain < matlab.mixin.CustomDisplay
 				
 				obj.grid = myParser.Results.grid(:);
 				obj.name = myParser.Results.name;
+				myParser.unmatchedWarning();
 			else
 				obj = quantity.Domain.empty();
 			end
-			
-			myParser.unmatchedWarning();
 		end
 		
 		function i = get.ones(obj)
@@ -82,45 +81,58 @@ classdef Domain < matlab.mixin.CustomDisplay
 			s = [obj.n];
 		end
 		
-		function s = eq(obj, obj2)
-			
-			if isempty(obj) && isempty(obj2)
-				% if both are empty -> they are equal
-				s = true;
-			elseif isempty(obj) || isempty(obj2)
-				% if only one domain is empty -> they are not equal
-				s = false;
-			else
-				% if both are not empty: check if the grids and the
-				% gridNames are equal
-				s = all(obj.gridLength == obj2.gridLength);
-			
-			for k = 1:numel(obj)	
-				s = s && all(obj(k).grid == obj2(k).grid);
-				s = s && strcmp(obj(k).name, obj2(k).name);
-			end
-			end
-		end
-		
+% 		function s = eq(obj, obj2)
+% 			
+% 			if isempty(obj) && isempty(obj2)
+% 				% if both are empty -> they are equal
+% 				s = true;
+% 			elseif isempty(obj) || isempty(obj2)
+% 				% if only one domain is empty -> they are not equal
+% 				s = false;
+% 			else
+% 				% if both are not empty: check if the grids and the
+% 				% gridNames are equal
+% 				s = all(obj.gridLength == obj2.gridLength);
+% 			
+% 			for k = 1:numel(obj)	
+% 				s = s && all(obj(k).grid == obj2(k).grid);
+% 				s = s && strcmp(obj(k).name, obj2(k).name);
+% 			end
+% 			end
+% 		end
+% 		
 		function l = ne(obj, obj2)
 			l = ~(obj == obj2);
 		end
 				
-		function [joinedDomain, index] = gridJoin(obj1, obj2)
+		function [joinedDomain, index] = join(domain1, domain2)
 			%% 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.
+			arguments
+				domain1 quantity.Domain = quantity.Domain.empty
+				domain2 quantity.Domain = quantity.Domain.empty
+			end
+			
+			if nargin == 1 && length(domain1) == 1
+				joinedDomain = domain1;
+				index = 1;
+				return;
+			end
 			
-			[joinedGrid, index] = unique({obj1.name, obj2.name}, 'stable');
+			[joinedGrid, index] = unique({domain1.name, domain2.name}, 'stable');
 			
 			joinedDomain = quantity.Domain.empty;
 			
+			% #todo@domain: the gird comparison seems to be very
+			% complicated
+			
 			% check for each grid if it is in the domain of obj1 or obj2 or
 			% both
 			for i = 1 : numel(joinedGrid)
 				currentGridName = joinedGrid{i};
-				[index1, logicalIndex1] = obj1.gridIndex(currentGridName);
-				[index2, logicalIndex2] = obj2.gridIndex(currentGridName);
+				[index1, logicalIndex1] = domain1.gridIndex(currentGridName);
+				[index2, logicalIndex2] = domain2.gridIndex(currentGridName);
 				
 				%
 				% 				if ~any(logicalIndex1)
@@ -132,8 +144,8 @@ classdef Domain < matlab.mixin.CustomDisplay
 				% Check if a domain is in both domains:
 				%	-> then take the finer one of both
 				if any(logicalIndex1) && any(logicalIndex2)
-					tempDomain1 = obj1(index1);
-					tempDomain2 = obj2(index2);
+					tempDomain1 = domain1(index1);
+					tempDomain2 = domain2(index2);
 					
 					assert(tempDomain1.lower == tempDomain2.lower, 'Grids must have same domain for gridJoin')
 					assert(tempDomain1.upper == tempDomain2.upper, 'Grids must have same domain for gridJoin')
@@ -145,17 +157,21 @@ classdef Domain < matlab.mixin.CustomDisplay
 					end
 				% If it is not in both, -> just take the normal grid
 				elseif any(logicalIndex1)
-					joinedDomain(i) = obj1(index1);
+					joinedDomain(i) = domain1(index1);
 				elseif any(logicalIndex2)
-					joinedDomain(i) = obj2(index2);
+					joinedDomain(i) = domain2(index2);
 				end
 				
 			end
+		end
+		
+		function [joinedDomain, index] = gridJoin(obj1, obj2)
+			[joinedDomain, index] = obj1.join(obj2);
 		end % gridJoin()
 				
 		function [idx, logicalIdx] = gridIndex(obj, names)
 			%% GRIDINDEX returns the index of the grid
-			% [idx, log] = gridIndex(obj, names) searches in the gridName
+			% [idx, log] = gridIndex(obj, names) searches in the name
 			% properties of obj for the "names" and returns its index as
 			% "idx" and its logical index as "log"
 			
@@ -174,7 +190,11 @@ classdef Domain < matlab.mixin.CustomDisplay
 				log = strcmp({obj.name}, names{k});
 				logicalIdx = logicalIdx | log;
 				if any(log)
-					idx(k) = nArgIdx(log);
+					% #todo@domain: if the index for a grid name is searched in a
+					% list with multiple grids but same names, only the first
+					% occurrence is returned. this should be changed...
+					tmpFirst = nArgIdx(log);
+					idx(k) = tmpFirst(1);
 				else
 					idx(k) = 0;
 				end
diff --git a/+quantity/Function.m b/+quantity/Function.m
index ac7157dafb72d8e80d563b897238c2534123bc37..dc15eb148e7cbe49544e23c5928fbe061294f343 100644
--- a/+quantity/Function.m
+++ b/+quantity/Function.m
@@ -81,7 +81,7 @@ classdef Function < quantity.Discrete
 			if nargin < 3
 				recalculate = false;
 			end
-			if ~recalculate && (myDomain == obj(1).domain)
+			if ~recalculate && isequal(myDomain, obj(1).domain)
 				value = obj2value@quantity.Discrete(obj);
 			else
 				% otherwise the function has to be evaluated on the new
diff --git a/+unittests/+quantity/testDomain.m b/+unittests/+quantity/testDomain.m
index c6b9738a8d0ad827eac392467d6054cd9c6ec302..5d899afef50fc2d6344b0532cc975489479950a8 100644
--- a/+unittests/+quantity/testDomain.m
+++ b/+unittests/+quantity/testDomain.m
@@ -5,15 +5,18 @@ end
 function setupOnce(testCase)
 end
 
-function testEq(testCase)
+function testEquality(testCase)
 
 d = quantity.Domain('grid', 0, 'name', 'd');
 e = quantity.Domain.empty();
 
-testCase.verifyFalse( d == e );
-testCase.verifyTrue( d == d );
-testCase.verifyTrue( e == e );
-testCase.verifyFalse( e == d );
+dd = d;
+ee = e;
+
+testCase.verifyNotEqual( d, e );
+testCase.verifyEqual( dd, d );
+testCase.verifyEqual( ee, e );
+testCase.verifyNotEqual( e, d );
 
 end