diff --git a/+misc/diagNd.m b/+misc/diagNd.m
index 88a3e555c84c24339ecb8d3d78177b5d1a21b725..263b7cb50da62503432749a579e9259b58e03452 100644
--- a/+misc/diagNd.m
+++ b/+misc/diagNd.m
@@ -23,50 +23,30 @@ function result = diagNd(value, dimensionsToBeDiagonalized)
 		dimensionsToBeDiagonalized = 1:ndims(value);
 	end % arguments
 	
-% 	sizeValue = size(value);
-% 	sizeResult = [sizeValue(dimensionsToBeDiagonalized(1)), ...
-% 		sizeValue(all((1:numel(sizeValue)) ~= dimensionsToBeDiagonalized.'))];
-% 	result = zeros(sizeResult);
-% 	
-% 	tempValueIdx = cell(numel(sizeValue), 1);
-% 	for it = 1 : numel(tempValueIdx)
-% 		if any(it == dimensionsToBeDiagonalized)
-% 			tempValueIdx{it} = 1;
-% 		else
-% 			tempValueIdx{it} = (1:1:sizeValue(it)).';
-% 		end
-% 	end
-% 	
-% 	for it = 1 : sizeResult(1)
-% 		for jt = 1 : numel(dimensionsToBeDiagonalized)
-% 			tempValueIdx{dimensionsToBeDiagonalized(jt)} = it;
-% 		end
-% 		tempValue = value(tempValueIdx{:});
-% 		result(it, :) = tempValue(:);
-% 	end
-	
-	% permute value such that all dimensions that are considered for
-	% diagonalization are at the beginning
-	originalSize = size(value);
-	originalOrder = 1:numel(originalSize);
-	value = permute(value, [dimensionsToBeDiagonalized, ...
-		originalOrder(all(originalOrder ~= dimensionsToBeDiagonalized(:)))]);
-
-	% the result contains the diagonalized dimensions first, then the
-	% remaining
-	tempValueSize = size(value);
-	sizeResult = [tempValueSize(numel(dimensionsToBeDiagonalized) : end), 1];
+	% save sizes and init result
+	sizeValue = size(value);
+	sizeResult = [sizeValue(dimensionsToBeDiagonalized(1)), ...
+		sizeValue(all((1:numel(sizeValue)) ~= dimensionsToBeDiagonalized.')), 1];
 	result = zeros(sizeResult);
 	
-% 	an n-dimensional identityMatrix = valueSelector which ensures 
-%		valueSelector(it, jt, kt, ...) = 1 if it = jt = kt = ...
-%									   = 0 else
-%	is used to select the diagonal elements in the first dimensions
-	valueSelector = misc.eyeNd(originalSize(dimensionsToBeDiagonalized), "logical");
+	% the cell array tempValueIdx will be used to read the elements of value.
+	tempValueIdx = cell(numel(sizeValue), 1);
+	for it = 1 : numel(tempValueIdx)
+		if any(it == dimensionsToBeDiagonalized)
+			tempValueIdx{it} = 1;
+		else
+			tempValueIdx{it} = (1:1:sizeValue(it)).';
+		end
+	end
 	
-	% collapse the first n dimensions to be diagonalized, so that the diagonal
-	% elements can be selected with the valueSelector by logical indexing efficiently
-	% afterwards.
-	valueCollapsed = reshape(value, [sizeResult(1)^(numel(dimensionsToBeDiagonalized)), sizeResult(2:end)]);
-	result(:,:) = valueCollapsed(valueSelector(:), :);
-end
+	% iterate over diagonal elements
+	for it = 1 : sizeResult(1)
+		% set same idx for all elements on diagonal
+		for jt = 1 : numel(dimensionsToBeDiagonalized)
+			tempValueIdx{dimensionsToBeDiagonalized(jt)} = it;
+		end
+		% read values from value and set values to result
+		tempValue = value(tempValueIdx{:});
+		result(it, :) = tempValue(:);
+	end
+end % diagNd()
\ No newline at end of file
diff --git a/+quantity/Discrete.m b/+quantity/Discrete.m
index e497a197d3d26ffe63b1258d9d020914ea685910..42e710e13833256f51a1b7f4072457fc9823398c 100644
--- a/+quantity/Discrete.m
+++ b/+quantity/Discrete.m
@@ -924,58 +924,57 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete ...
 						% f(zeta, zeta) -> the 2nd subs(f, zeta, z) will
 						% result in f(z, z) and not in f(zeta, z) as
 						% intended. This is solved, by an additonal
-						% substitution:
+						% substitution: 
+						%	f.subs(z,zetabackUp).subs(zeta,z).subs(zetabackUp,zeta)
 						values{end+1} = values{1};
 						gridName2Replace{end+1} = [gridName2Replace{1}, 'backUp'];
-						values{1} = [gridName2Replace{1}, 'backUp'];
+						values{1} = gridName2Replace{end};
 					end
 					if isequal(values{1}, gridName2Replace{1})
-						% replace with same variable... everything stay the
+						% replace with same variable... everything stays the
 						% same.
-						newGrid = obj(1).grid;
-						newGridName = obj(1).gridName;
+						% Do not use "return", since, later subs might need to be
+						% called recursively!
 						newValue = obj.on();
+						newDomain = obj(1).domain;
 					elseif any(strcmp(values{1}, obj(1).gridName))
 						% if for a quantity f(z, zeta) this method is
 						% called with subs(f, zeta, z), then g(z) = f(z, z)
 						% results, hence the dimensions z and zeta are
 						% merged.
-						gridIndices = [obj(1).domain.gridIndex(gridName2Replace{1}), ...
+						domainIndices = [obj(1).domain.gridIndex(gridName2Replace{1}), ...
 							obj(1).domain.gridIndex(values{1})];
-						newGridForOn = obj(1).grid;
-						if numel(obj(1).grid{gridIndices(1)}) > numel(obj(1).grid{gridIndices(2)})
-							newGridForOn{gridIndices(2)} = newGridForOn{gridIndices(1)};
-						else
-							newGridForOn{gridIndices(1)} = newGridForOn{gridIndices(2)};
+						newDomainForOn = obj(1).domain;
+						if obj(1).domain(domainIndices(1)).n > obj(1).domain(domainIndices(2)).n
+							newDomainForOn(domainIndices(2)) = quantity.Domain(...
+								newDomainForOn(domainIndices(2)).name, ...
+								newDomainForOn(domainIndices(1)).grid);
+						elseif  obj(1).domain(domainIndices(1)).n < obj(1).domain(domainIndices(2)).n
+							newDomainForOn(domainIndices(1)) = quantity.Domain(...
+								newDomainForOn(domainIndices(1)).name, ...
+								newDomainForOn(domainIndices(2)).grid);
 						end
-						newValue = misc.diagNd(obj.on(newGridForOn), gridIndices);
-						newGrid = {newGridForOn{gridIndices(1)}, ...
-							newGridForOn{1:1:numel(newGridForOn) ~= gridIndices(1) ...
-							& 1:1:numel(newGridForOn) ~= gridIndices(2)}};
-						newGridName = {values{1}, ...
-							obj(1).gridName{1:1:numel(obj(1).gridName) ~= gridIndices(1) ...
-							& 1:1:numel(obj(1).gridName) ~= gridIndices(2)}};
+						newValue = misc.diagNd(obj.on(newDomainForOn), domainIndices);
+						newDomain = [newDomainForOn(domainIndices(2)), ...
+							newDomainForOn(all(1:1:numel(newDomainForOn) ~= domainIndices(:)))];
 					else
-						% this is the default case. just grid name is
-						% changed.
-						newGrid = obj(1).grid;
-						newGridName = obj(1).gridName;
-						newGridName{obj(1).domain.gridIndex(gridName2Replace{1})} ...
-							= values{1};
+						% this is the default case. just grid name is changed.
+						newDomain = obj(1).domain;
+						newDomain(obj(1).domain.gridIndex(gridName2Replace{1})) = ...
+							quantity.Domain(values{1}, ...
+							obj(1).domain(obj(1).domain.gridIndex(gridName2Replace{1})).grid);
 						newValue = obj.on();
 					end
 					
 				elseif isnumeric(values{1}) && numel(values{1}) == 1
 					% if values{1} is a scalar, then obj is evaluated and
-					% the resulting quantity loses that spatial grid and
+					% the resulting quantity looses that spatial grid and
 					% gridName
-					newGridName = obj(1).gridName;
-					newGridName = newGridName(~strcmp(newGridName, gridName2Replace{1}));
+					newDomain = obj(1).domain;
+					newDomain = newDomain(~strcmp({newDomain.name}, gridName2Replace{1}));
 					% newGrid is the similar to the original grid, but the
 					% grid of gridName2Replace is removed.
-					newGrid = obj(1).grid;
-					newGrid = newGrid((1:1:numel(newGrid)) ~= obj(1).domain.gridIndex(gridName2Replace{1}));
-					newGridSize = cellfun(@(v) numel(v), newGrid);
+					newGridSize = newDomain.gridLength();
 					% newGridForOn is the similar to the original grid, but
 					% the grid of gridName2Replace is set to values{1} for
 					% evaluation of obj.on().
@@ -986,18 +985,18 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete ...
 				elseif isnumeric(values{1}) && numel(values{1}) > 1
 					% if values{1} is a double vector, then the grid is
 					% replaced.
-					newGrid = obj(1).grid;
-					newGrid{obj(1).domain.gridIndex(gridName2Replace{1})} = values{1};
-					newGridName = obj(1).gridName;
-					newValue = obj.on(newGrid);
+					newDomain = obj(1).domain;
+					newDomain(obj(1).domain.gridIndex(gridName2Replace{1})) = ...
+						quantity.Domain(gridName2Replace{1}, values{1});
+					newValue = obj.on(newDomain);
 				else
 					error('value must specify a gridName or a gridPoint');
 				end
-				if isempty(newGridName)
+				if isempty(newDomain)
 					solution = newValue;
 				else
 					solution = quantity.Discrete(newValue, ...
-						'grid', newGrid, 'gridName', newGridName, ...
+						'domain', newDomain, ...
 						'name', obj(1).name);
 				end
 				if numel(gridName2Replace) > 1
@@ -2142,7 +2141,7 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete ...
 			
 			value = reshape(v, [obj(1).domain.gridLength(), size(obj)]);
 				
-			if nargin == 2 && ~isequal( 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
diff --git a/+quantity/Domain.m b/+quantity/Domain.m
index e58be9d42912b0d800152a9a60dc3ef30b29fb66..a9669ef4ca904c500de060f10bef1705f385674c 100644
--- a/+quantity/Domain.m
+++ b/+quantity/Domain.m
@@ -7,13 +7,13 @@ classdef (InferiorClasses = {?quantity.EquidistantDomain}) Domain < ...
 		% continuous quantity. For an example, the function f(x) should be
 		% considered on the domain x \in X = [0, 1]. Then, a grid can be
 		% generated by X_grid = linspace(0, 1).
-		grid double {mustBeReal, mustBeFinite};
+		grid (:, 1) double;
 		
 		% a speaking name for this domain; Should be unique, so that the
 		% domain can be identified by the name.
 		name string = "";
 		
-		n (1,1) double {mustBeInteger};		% number of discretization points == gridLength
+		n (1,1) double;		% number of discretization points == gridLength
 	end
 	
 	properties (Dependent)
@@ -30,7 +30,7 @@ classdef (InferiorClasses = {?quantity.EquidistantDomain}) Domain < ...
 			%DOMAIN initialize the domain
 			%
 			if nargin >= 1
-				obj.grid = grid(:);
+				obj.grid = grid;
 				obj.name = name;
 				obj.n = length(obj.grid);
 			end