diff --git a/+misc/ensureString.m b/+misc/ensureString.m
new file mode 100644
index 0000000000000000000000000000000000000000..52d80b02d7e59fb9a7cfcc75e64df49265c32bd9
--- /dev/null
+++ b/+misc/ensureString.m
@@ -0,0 +1,16 @@
+function [str] = ensureString(chr)
+
+if iscell( chr )
+	if ~all( cellfun( @ischar, chr ) )
+		chr = [chr{:}];
+	end
+end
+
+if isstring( chr )
+	str = chr;
+else
+	str = convertCharsToStrings( chr );
+end
+
+end
+
diff --git a/+mustBe/gridName.m b/+mustBe/gridName.m
new file mode 100644
index 0000000000000000000000000000000000000000..d2c4147894615138cca664783c65cb21fb1f3bbb
--- /dev/null
+++ b/+mustBe/gridName.m
@@ -0,0 +1,8 @@
+function gridName( name )
+
+name = misc.ensureString( name );
+
+assert( isstring( name ), 'A grid name should be a string-array')
+
+end
+
diff --git a/+quantity/Discrete.m b/+quantity/Discrete.m
index b435a3cfb2daeb688ed1b1297839789047f55ebc..63841fcfd54642f20eacb2e04ff528934e930c64 100644
--- a/+quantity/Discrete.m
+++ b/+quantity/Discrete.m
@@ -1,4 +1,5 @@
-classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mixin.Copyable & matlab.mixin.CustomDisplay
+classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete ...
+		< handle & matlab.mixin.Copyable & matlab.mixin.CustomDisplay
 
 	properties (SetAccess = protected)
 		% Discrete evaluation of the continuous quantity
@@ -79,7 +80,7 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
 				
 				%% input parser
 				myParser = misc.Parser();
-				myParser.addParameter('name', string(), @isstr);
+				myParser.addParameter('name', "", @mustBe.gridName);
 				myParser.addParameter('figureID', 1, @isnumeric);
 				myParser.parse(varargin{:});
 				
@@ -143,9 +144,9 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
 		%---------------------------
 		function gridName = get.gridName(obj)
 			if isempty(obj.domain)
-				gridName = {};
+				gridName = [];
 			else
-				gridName = {obj.domain.name};
+				gridName = [obj.domain.name];
 			end
 		end
 		
@@ -316,21 +317,18 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
 			%	for this, find the index of the common domain in list of
 			%	temporary combined domain
 			
-			intersectDomain = intersect( {originalDomain( ~logOfDomain ).name}, ...
-				{g(1).domain.name} );
+			intersectDomain = intersect( originalDomain( ~logOfDomain ), ...
+				g(1).domain );
 			
 			if ~isempty(intersectDomain)
 				
-				idx = 1:length(tmpDomain);
-				logCommon = strcmp({tmpDomain.name}, intersectDomain);
+				idx = tmpDomain.gridIndex( intersectDomain );
 				
 				% take the diagonal values of the common domain, i.e., z = zeta		
 				% use the diag_nd function because it seems to be faster
 				% then the diagNd function, although the values must be
 				% sorted.
-				% #TODO: Rewrite the diagNd function, using for loops, in order to be as fast as diag_nd
-				newValues = permute( newValues, [idx(logCommon), idx(~logCommon)]);
-				newValues = misc.diag_nd(newValues);
+				newValues = misc.diagNd(newValues, idx);
 			end
 			
 			% *) build a new valueDiscrete on the correct grid.		
@@ -398,14 +396,10 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
 						% Since the order of the domains is not neccessarily equal to the
 						% order in obj(1).domain, this is more involved:
 						myDomain = misc.ensureIsCell(myDomain);
-						gridNames = misc.ensureIsCell(gridNames);
+						gridNames = misc.ensureString(gridNames);
 
 						assert(all(cellfun(@(v)isvector(v), myDomain)), ...
 							'The cell entries for a new grid have to be vectors')
-						assert(iscell(gridNames), ...
-							'The gridNames parameter must be cell array')
-						assert(all(cellfun(@ischar, gridNames)), ...
-							'The gridNames must be strings')
 
 						newGrid = myDomain;
 						myDomain = quantity.Domain.empty();
@@ -1015,7 +1009,7 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
 		
 		function [idx, logicalIdx] = gridIndex(obj, varargin)
 			warning('DEPRICATED: use quantity.Domain.gridIndex method instead')
-			[idx, logicalIdx] = obj(1).domain.gridIndex(varargin{:});
+			[~, idx, logicalIdx] = obj(1).domain.find(varargin{:});
 		end 
 		
 		function value = at(obj, point)
@@ -1305,19 +1299,20 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
 			%	corresponding domain from DOMAIN.
 			%
 			%	example:
-			%	q.changeGrid( linspace(0,1)', 't')
-			%	will change the grid with the name 't' to the new grid linspace(0,1)'
+			%		q.changeGrid( linspace(0,1)', 't')
+			%	will change the grid with the name 't' to the new grid
+			%	linspace(0,1)'
 			if isempty(obj)
 				newObj = obj.copy();
 				return;
 			end
 			
 			if isa(gridNew, 'quantity.Domain')
-				gridNameNew = {gridNew.name};
+				gridNameNew = [gridNew.name];
 				gridNew = {gridNew.grid};
 			else
+				gridNameNew = misc.ensureString(gridNameNew);				
 				gridNew  = misc.ensureIsCell(gridNew);
-				gridNameNew = misc.ensureIsCell(gridNameNew);
 			end
 			
 			if obj(1).isConstant
@@ -1325,7 +1320,7 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
 				newDomain(1:length( gridNew )) = quantity.Domain();
 				for it = 1 : length(gridNew)
 					newDomain(it) = ...
-						quantity.Domain(gridNameNew{it}, gridNew{it});
+						quantity.Domain(gridNameNew(it), gridNew{it});
 				end
 			else
 				gridIndexNew = obj(1).domain.gridIndex(gridNameNew);
@@ -1333,9 +1328,9 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
 
 				for it = 1 : length(gridIndexNew)
 					newDomain(gridIndexNew(it)) = ...
-						quantity.Domain(gridNameNew{it}, gridNew{it});
+						quantity.Domain(gridNameNew(it), gridNew{it});
 				end
-				assert(isequal({newDomain.name}, obj(1).gridName), ...
+				assert(isequal([newDomain.name], obj(1).gridName), ...
 					'rearranging grids failed');
 			end
 			
@@ -1661,17 +1656,13 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
 			% Optionally, a weight can be defined, such that instead
 			%	xNorm = sqrt(int_0^1 x.' * weight * x dz).
 			% The integral domain is specified by integralGrid.
-			
 			arguments
 				obj;
-				integralGridName = {obj(1).domain.name};
+				integralGridName {mustBe.gridName} = [obj(1).domain.name];
 				optArg.weight = eye(size(obj, 1));
 			end
 			
-			integralGridName = misc.ensureIsCell(integralGridName);
-			
-			assert(all(ischar([integralGridName{:}])), ...
-				'integralGrid must specify a gridName as a char');
+			integralGridName = misc.ensureString(integralGridName);
 			
 			if obj.nargin == 1 && all(strcmp(obj(1).gridName, integralGridName))
 				xNorm = sqrtm(on(int(obj.' * optArg.weight * obj), integralGridName));
@@ -1844,16 +1835,16 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
 			end
 		end
 		
-		function result = diff(obj, k, diffGridName)
+		function result = diff(obj, diffGridName, k)
 			% DIFF computation of the derivative
 			%	result = DIFF(obj, k, diffGridName) applies the
 			% 'k'th-derivative for the variable specified with the input
 			% 'diffGridName' to the obj. If no 'diffGridName' is specified,
 			% then diff applies the derivative w.r.t. all gridNames.
-			if nargin == 1 || isempty(k)
-				k = 1;  % by default, only one derivatve per diffGridName is applied
-			else
-				assert(isnumeric(k) && (round(k) == k))
+			arguments
+				obj;
+				diffGridName (1,1) string = obj(1).gridName;
+				k uint64 = 1;
 			end
 			
 			if obj.isConstant && isempty(obj(1).gridName)
@@ -1864,11 +1855,6 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
 				return
 			end
 			
-			if nargin < 3 % if no diffGridName is specified, then the derivative
-				% w.r.t. all gridNames is applied
-				diffGridName = obj(1).gridName;
-			end
-			
 			if isa(diffGridName, 'quantity.Domain')
 				% a quantity.Domain is used instead of a grid name 
 				%	-> get the grid name from the domain object
@@ -1880,22 +1866,7 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
 				diffGridName = {diffGridName.name};
 			end
 			
-			% diff for each element of diffGridName (this is rather
-			% inefficient, but an easy implementation of the specification)
-			if iscell(diffGridName) || isempty(diffGridName)
-				if numel(diffGridName) == 0 || isempty(diffGridName)
-					result = copy(obj);
-				else
-					result = obj.diff(k, diffGridName{1}); % init result
-					for it = 2 : numel(diffGridName)
-						result = result.diff(k, diffGridName{it});
-					end
-				end
-			elseif k == 0
-				result = obj;
-			else
-				result = obj.diff_inner(k, diffGridName);
-			end
+			result = obj.diff_inner(k, diffGridName);
 		end
 		
 		function I = int(obj, varargin)
@@ -2323,8 +2294,8 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
 		end
 		
 		function result = diff_inner(obj, k, diffGridName)
-			gridSelector = strcmp(obj(1).gridName, diffGridName);
-			gridSelectionIndex = find(gridSelector);
+			
+			gridSelectionIndex = obj(1).domain.gridIndex(diffGridName);
 
 			spacing = gradient(obj(1).grid{gridSelectionIndex}, 1);
 			assert(numeric.near(spacing, spacing(1)), ...
@@ -2350,9 +2321,9 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
 				'name', ['(d_{', diffGridName, '}', obj(1).name, ')']);
 
 			if k > 1
-				% 			% if a higher order derivative is requested, call the function
-				% 			% recursivly until the first-order derivative is reached
-				result = result.diff(k-1, diffGridName);
+				% if a higher order derivative is requested, call the function
+				% recursivly until the first-order derivative is reached
+				result = result.diff(diffGridName, k-1);
 			end
 		end
 		
@@ -2360,13 +2331,18 @@ classdef  (InferiorClasses = {?quantity.Symbolic}) Discrete < handle & matlab.mi
 			% Computes the required permutation vectors to use
 			% misc.multArray for multiplication of matrices
 			
-			uA = unique(a(1).gridName, 'stable');
-			assert(numel(uA) == numel(a(1).gridName), 'Gridnames have to be unique!');
-			uB = unique(b(1).gridName, 'stable');
-			assert(numel(uB) == numel(b(1).gridName), 'Gridnames have to be unique!');
+			% #todo this assertion should not be required!
+% 			uA = unique(a(1).gridName, 'stable');
+% 			assert(numel(uA) == numel(a(1).gridName), 'Gridnames have to be unique!');
+% 			uB = unique(b(1).gridName, 'stable');
+% 			assert(numel(uB) == numel(b(1).gridName), 'Gridnames have to be unique!');
 			
 			% 1) find common entries
-			common = intersect(a(1).gridName, b(1).gridName);
+			if isempty( b(1).gridName ) || isempty( a(1).gridName )
+				common = [];
+			else
+				common = intersect(a(1).gridName, b(1).gridName);
+			end
 			
 			commonA = false(1, a.nargin());
 			commonB = false(1, b.nargin());
diff --git a/+quantity/Domain.m b/+quantity/Domain.m
index 57170254dc441170ae05520589727c877fe5a66d..1e10a1aaacb5f8ffd67d1841a5110fb7910475c9 100644
--- a/+quantity/Domain.m
+++ b/+quantity/Domain.m
@@ -11,7 +11,7 @@ classdef (InferiorClasses = {?quantity.EquidistantDomain}) Domain < ...
 		
 		% a speaking name for this domain; Should be unique, so that the
 		% domain can be identified by the name.
-		name char;
+		name string = "";
 		
 		n (1,1) double {mustBeInteger};		% number of discretization points == gridLength
 	end
@@ -96,7 +96,7 @@ classdef (InferiorClasses = {?quantity.EquidistantDomain}) Domain < ...
 			s = [obj.n];
 		end
 		
-		function d = find(obj, domain, varargin)
+		function d = find(obj, searchName)
 			% FIND a domain in object-array of quantity.Domain
 			%
 			%	d = find(OBJ, NAME) returns the domain with the name NAME
@@ -111,30 +111,56 @@ classdef (InferiorClasses = {?quantity.EquidistantDomain}) Domain < ...
 			%
 			%	d = find(OBJ, { NAME or DOMAIN}) the NAME or the DOMAIN to
 			%	be found can also be defined by cell-array.
-			
-			names = {obj.name};
-			if isa(domain, 'quantity.Domain')
-				searchName = {domain.name};
-			elseif ischar(domain)
-				searchName = {domain};
-			elseif iscell(domain)
-				searchName = domain;
-			else
-				error('Not implemented yet')
+			arguments
+				obj quantity.Domain;
+			end
+			arguments (Repeating)
+				searchName string;
+			end
+						
+			d = obj( obj.gridIndex( searchName ) );
+		end
+		function [idx, log] = gridIndex(obj, searchName, position)
+			%% GRIDINDEX returns the index of the grid
+			% [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"
+			arguments
+				obj
+				searchName = [obj.name];
+				position string = "all";
 			end
 			
-			if numel(searchName) > 1
-				d = obj.find(searchName{:});
-				return
+			if isa(searchName, 'quantity.Domain')
+				searchName = [searchName.name];
 			end
 			
-			d = obj(strcmp(names, searchName));
+			objNames = [obj.name];
+			
+			if isempty(objNames)
+				idx = 0;
+				log = [];
+			else
+				searchName = misc.ensureString( searchName );
+				
+				log = false(size(objNames));
+				counter = 1:numel(objNames);
+				idx = [];
+				for k = 1 : numel(searchName)
+					log_ = objNames == searchName(k);
+					log = log | log_;
+					idx = [idx counter( log_ )];
+				end
+				if isempty(idx)
+					idx = 0;
+				end
+			end
 			
-			for it = 1 : numel(varargin)
-				d = [d, find(obj, varargin{:})];
+			if position == "first"
+				idx = idx(1);
 			end
-		end
-		
+			
+		end % gridIndex		
 		function l = ne(obj, obj2)
 			% ~= Not equal.
 			l = ~(obj == obj2);
@@ -155,27 +181,17 @@ classdef (InferiorClasses = {?quantity.EquidistantDomain}) Domain < ...
 				return;
 			end
 			
-			[joinedGrid, index] = unique({domain1.name, domain2.name}, 'stable');
+			[joinedGrid, index] = unique([domain1.name, domain2.name], 'stable');
 			
 			joinedDomain(1 : numel(joinedGrid)) = quantity.Domain();
 			
-			% #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] = domain1.gridIndex(currentGridName);
-				[index2, logicalIndex2] = domain2.gridIndex(currentGridName);
-				
-				%
-				% 				if ~any(logicalIndex1)
-				% 					joinedDomain(i) = obj2(index2);
-				% 				elseif ~any(logicalIndex2)
-				% 					joinedDomain(i) = obj1(index1);
-				% 				else
-				%
+				currentGridName = joinedGrid(i);
+				[index1, logicalIndex1] = domain1.gridIndex(currentGridName, "first");
+				[index2, logicalIndex2] = domain2.gridIndex(currentGridName, "first");
+								
 				% Check if a domain is in both domains:
 				%	-> then take the finer one of both
 				if any(logicalIndex1) && any(logicalIndex2)
@@ -192,7 +208,7 @@ classdef (InferiorClasses = {?quantity.EquidistantDomain}) Domain < ...
 					end
 				% If it is not in both, -> just take the normal grid
 				elseif any(logicalIndex1)
-					joinedDomain(i) = domain1(index1);
+					joinedDomain(i) = domain1(index1(1));
 				elseif any(logicalIndex2)
 					joinedDomain(i) = domain2(index2);
 				end
@@ -204,41 +220,7 @@ classdef (InferiorClasses = {?quantity.EquidistantDomain}) Domain < ...
 			[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 name
-			% properties of obj for the "names" and returns its index as
-			% "idx" and its logical index as "log"
-			arguments
-				obj
-				names = {obj.name};
-			end
-			
-			if isa(names, 'quantity.Domain')
-				names = {names.name};
-			end
-			
-			names = misc.ensureIsCell(names);
-			
-			idx = zeros(1, length(names));
-			nArgIdx = 1:obj.ndims();
-			
-			logicalIdx = false(1, obj.ndims());
-			
-			for k = 1:length(names)
-				log = strcmp({obj.name}, names{k});
-				logicalIdx = logicalIdx | log;
-				if any(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
-			end
-		end % gridIndex
+
 		
 		function i = isempty(obj)
 			i = any(size(obj) == 0);
@@ -269,14 +251,14 @@ classdef (InferiorClasses = {?quantity.EquidistantDomain}) Domain < ...
 		function [idx, newDomain] = getPermutationIdx(obj, order)
 			
 			if isa(order, 'quantity.Domain')
-				names = {order.name};
-			elseif ischar([order{:}])
-				names = order;
+				names = [order.name];
+			elseif iscell(order) || ischar(order) || isstring(order)
+				names = misc.ensureString( order );
 			else
-				error('the input parameter order must be a array of quantity.Domain objects or a cell-array with string')
+				error('the input parameter order must be a array of quantity.Domain objects or a string-array')
 			end
 						
-			idx = cellfun(@(v) obj.gridIndex(v), names);
+			idx = cellfun(@(v) obj.gridIndex(v), names); % #todo@domainNameString
 			
 			if isa(order, 'quantity.Domain')
 				newDomain = order;
@@ -302,7 +284,7 @@ classdef (InferiorClasses = {?quantity.EquidistantDomain}) Domain < ...
 			
 			% only sort the grids if there is something to sort
 			if obj.ndims > 1
-				gridNames = {obj.name};
+				gridNames = [obj.name];
 				
 				% this is the default case for ascending alphabetical
 				% order
@@ -404,7 +386,7 @@ classdef (InferiorClasses = {?quantity.EquidistantDomain}) Domain < ...
 			%% domain parser
 			domainParser = misc.Parser();
 			domainParser.addParameter('domain', {}, @(g) isa(g, 'quantity.Domain'));
-			domainParser.addParameter('gridName', '', @(g) ischar(g) || iscell(g));
+			domainParser.addParameter('gridName', '');
 			domainParser.addParameter('grid', [], @(g) isnumeric(g) || iscell(g));
 			domainParser.parse(varargin{:});
 
@@ -418,7 +400,7 @@ classdef (InferiorClasses = {?quantity.EquidistantDomain}) Domain < ...
 				%	-> initialize quantity.Domain objects with the
 				%	specified values
 
-				myGridName = misc.ensureIsCell(domainParser.Results.gridName);
+				myGridName = misc.ensureString(domainParser.Results.gridName);
 				myGrid = misc.ensureIsCell(domainParser.Results.grid);
 
 				assert(isequal(numel(myGrid), numel(myGridName)), ...
@@ -427,7 +409,7 @@ classdef (InferiorClasses = {?quantity.EquidistantDomain}) Domain < ...
 				% initialize the domain objects
 				myDomain = quantity.Domain.empty();
 				for k = 1:numel(myGrid)
-					myDomain(k) = quantity.Domain(myGridName{k}, myGrid{k});
+					myDomain(k) = quantity.Domain(myGridName(k), myGrid{k});
 				end
 			else
 				% else case: the domains are specified as domain
@@ -435,7 +417,7 @@ classdef (InferiorClasses = {?quantity.EquidistantDomain}) Domain < ...
 				myDomain = domainParser.Results.domain;
 			end	
 			
-			assert(numel(myDomain) == numel(unique({myDomain.name})), ...
+			assert(misc.isunique([myDomain.name]), ...
 				'The names of the domain must be unique');
 			
 			unmatched = domainParser.UnmatchedNameValuePair;
diff --git a/+quantity/Piecewise.m b/+quantity/Piecewise.m
index 4f2844d7c5772bc4ce7e18ba7dea1a44cc3d2d96..adfb6129fa4c01dcd6fb319fc6100855033bfc64 100644
--- a/+quantity/Piecewise.m
+++ b/+quantity/Piecewise.m
@@ -31,21 +31,21 @@ classdef (InferiorClasses = {?quantity.Symbolic}) Piecewise < quantity.Discrete
 
 				% ensure all have the same domain
 				d = quantities{1}(1).domain;
-				assert( all( cellfun(@(q) all( strcmp( {d.name}, {q(1).domain.name} ) ), quantities) ) , ...
+				assert( all( cellfun(@(q) all( strcmp( [d.name], [q(1).domain.name] ) ), quantities) ) , ...
 					'The quantities for the piecewise combination must have the same domain names' );
 				assert(length(quantities) >= 2, 'Only one quantity is given for the piecewise combination. At least 2 are required.');
 
 				assert( quantities{1}(1).domain.gridIndex(domainToJoin) == 1, ...
 					'The domain to join must be the first domain!');
 
-				joinedGrid = quantities{1}(1).domain.find(domainToJoin).grid;
+				joinedGrid = quantities{1}(1).domain.find(domainToJoin.name).grid;
 				joinedValues = quantities{1}.on();
 
 				% ensure the domains fit to each other and then concatenate
 				% them
 				for k = 2:length(quantities)
-					dkm1 = quantities{k-1}(1).domain.find(domainToJoin);
-					dk = quantities{k}(1).domain.find(domainToJoin);
+					dkm1 = quantities{k-1}(1).domain.find(domainToJoin.name);
+					dk = quantities{k}(1).domain.find(domainToJoin.name);
 
 					assert(dkm1.upper == dk.lower, 'Domains do not connect to each other');
 
diff --git a/+quantity/Symbolic.m b/+quantity/Symbolic.m
index db27771990d1fb852f45951737a8cfcd6ff1efe2..ae435a79832d76c80015975349639df3b2d9c338 100644
--- a/+quantity/Symbolic.m
+++ b/+quantity/Symbolic.m
@@ -70,10 +70,10 @@ classdef Symbolic < quantity.Function
 					% variable
 					if isa(fun{k}, 'sym')
 						symb{k} = fun{k};
-						fun{k} = quantity.Symbolic.setValueContinuous(symb{k}, {myDomain.name});
+						fun{k} = quantity.Symbolic.setValueContinuous(symb{k}, [myDomain.name]);
 					elseif isa(fun{k}, 'double')
 						symb{k} = sym(fun{k});
-						fun{k} = quantity.Symbolic.setValueContinuous(symb{k}, {myDomain.name});						
+						fun{k} = quantity.Symbolic.setValueContinuous(symb{k}, [myDomain.name]);						
 					elseif isa(fun{k}, 'function_handle')
 						symb{k} = sym(fun{k});
 					else
@@ -299,10 +299,10 @@ classdef Symbolic < quantity.Function
 				else
 					gridName = {gridName.name};
 				end
-			else
-				gridName = misc.ensureIsCell(gridName);
 			end
 			
+			gridName = misc.ensureString(gridName);
+			
 			if nargin == 3 && isa(values, 'quantity.Domain')
 				% replacement of the grid AND the gridName
 				% 1) replace the grid
@@ -314,14 +314,6 @@ classdef Symbolic < quantity.Function
 				values = misc.ensureIsCell(values);
 			end
 			
-% 			% ensure that domains which should be replaced are known by the
-% 			% object
-% 			for k = 1:length(values)
-% 				if ischar(values{k})
-% 					assert( any( strcmp( obj(1).gridName, values{k} ) ), 'The domain name to be replaced must be a domain name of the object' );
-% 				end
-% 			end
-			
 			isNumericValue = cellfun(@isnumeric, values);
 			if any((cellfun(@(v) numel(v(:)), values)>1) & isNumericValue)
 				error('only implemented for one value per grid');
@@ -344,13 +336,13 @@ classdef Symbolic < quantity.Function
  				for it = 1 : numel(values)
  					if ~isNumericValue(it) && ~isempty(obj(1).gridName(~selectRemainingGrid)) 
 						% check if there is a symbolic value and if this value exists in the object
-						if ischar(values{it})
-							newGridName{strcmp(obj(1).gridName(~selectRemainingGrid), gridName{it})} ...
+						if isstring( values{it} ) || ischar( values{it} )
+							newGridName( strcmp(obj(1).gridName(~selectRemainingGrid), gridName(it)) ) ...
 								= values{it};
 						elseif isa(values{it}, 'sym')
 							gridNameTemp = symvar(values{it});
 							assert(numel(gridNameTemp) == 1, 'replacing one gridName with 2 gridName is not supported');
-							newGridName{strcmp(obj(1).gridName(~selectRemainingGrid), gridName{it})} ...
+							newGridName{strcmp(obj(1).gridName(~selectRemainingGrid), gridName(it))} ...
 								= char(gridNameTemp);
 							
 						else
@@ -359,7 +351,7 @@ classdef Symbolic < quantity.Function
  					end
 				end
 								
-				symbolicSolution = subs(obj.sym(), gridName, values);				
+				symbolicSolution = subs(obj.sym(), cellstr(gridName), values);				
 				
 				if isempty(symvar(symbolicSolution(:)))
 					% take the new values that are not numeric as the new
@@ -374,7 +366,7 @@ classdef Symbolic < quantity.Function
 					end
 					for it = 1 : numel(newGridName)
 						if ~isempty(charValues) && any(contains(charValues, newGridName{it}))
-							newGrid{it} = obj.gridOf(gridName{strcmp(values, newGridName{it})});
+							newGrid{it} = obj.gridOf(gridName( strcmp(values, newGridName{it}) ));
 						else
 							newGrid{it} = obj.gridOf(newGridName{it});
 						end
@@ -709,11 +701,8 @@ classdef Symbolic < quantity.Function
 			%  from a to b."
 			
 			if nargin == 2
-				if ~iscell(z)
-					z = {z};
-				end
 				
-				grdIdx = strcmp(obj(1).gridName, z);
+				grdIdx = obj(1).domain.gridIndex(z);
 				a = obj(1).grid{grdIdx}(1);
 				b = obj(1).grid{grdIdx}(end);
 			end
@@ -949,7 +938,11 @@ classdef Symbolic < quantity.Function
 					fDouble = double(f);
 					f = @(varargin) fDouble + quantity.Function.zero(varargin{:});
 				else
-					f = matlabFunction(f, 'Vars', symVar );
+					if iscell( symVar )
+						symVar = [symVar{:}];
+					end
+						
+					f = matlabFunction(f, 'Vars', convertStringsToChars( symVar ) );
 				end
 			end
 		end % setValueContinuous()
diff --git a/+unittests/+misc/testMisc.m b/+unittests/+misc/testMisc.m
index 009f05ba264785a8fe5b371a2518cbffb87e3684..c8677b5826f9915e1eb74792ff4db25bf3219833 100644
--- a/+unittests/+misc/testMisc.m
+++ b/+unittests/+misc/testMisc.m
@@ -4,6 +4,18 @@ function [tests] = testMisc()
 tests = functiontests(localfunctions());
 end
 
+function testEnsureString(testCase)
+
+a = "a";
+b = 'b';
+
+testCase.verifyEqual( misc.ensureString(a), a)
+testCase.verifyEqual( misc.ensureString(b), string(b))
+testCase.verifyEqual( misc.ensureString({a, b}), [a, string(b)])
+testCase.verifyEqual( misc.ensureString({b, b}), [string(b), string(b)])
+
+end
+
 function testEyeNd(tc)
 n = 2;
 tc.verifyEqual(misc.eyeNd([n, n]), eye(n));
diff --git a/+unittests/+quantity/testDiscrete.m b/+unittests/+quantity/testDiscrete.m
index 2a45b978b129fc99a707a451eade0b8ecc733506..7cbe9f6959f2988c3a2cbd0776b1f5305939eb10 100644
--- a/+unittests/+quantity/testDiscrete.m
+++ b/+unittests/+quantity/testDiscrete.m
@@ -1,14 +1,27 @@
-
 function [tests] = testDiscrete()
 %testQuantity Summary of this function goes here
 %   Detailed explanation goes here
 tests = functiontests(localfunctions);
 end
 
+
+function setupOnce(testCase)
+syms z zeta t sy
+testCase.TestData.sym.z = z;
+testCase.TestData.sym.zeta = zeta;
+testCase.TestData.sym.t = t;
+testCase.TestData.sym.sy = sy;
+end
+
 function testL2Norm(tc)
+
+Z = tc.TestData.sym.z;
+T = tc.TestData.sym.t;
+ZETA = tc.TestData.sym.zeta;
+
 t = quantity.Domain('t', linspace(0, 2, 101));
 z = quantity.Domain('z', linspace(0, 1, 51));
-x = quantity.Symbolic([sym('z') * sym('t'); sym('t')], 'domain', [z, t]);
+x = quantity.Symbolic([Z * T; T], 'domain', [z, t]);
 
 xNorm = sqrt(int(x.' * x, 'z'));
 tc.verifyEqual(MAX(abs(xNorm - x.l2norm('z'))), 0, 'AbsTol', 1e-12);
@@ -105,9 +118,6 @@ tc.verifyEqual(q1V , q2V')
 end % testSort
 
 function testBlkdiag(tc)
-% init some data
-syms z zeta
-
 z = quantity.Domain('z', linspace(0, 1, 21));
 zeta = quantity.Domain('zeta', linspace(0, 1, 41));
 
@@ -149,7 +159,9 @@ tc.verifyEqual(blub(2).name, 'asdf');
 end
 
 function testCtranspose(tc)
-syms z zeta
+z = tc.TestData.sym.z;
+zeta = tc.TestData.sym.zeta;
+
 qSymbolic = quantity.Symbolic(...
 	[1+z*zeta, -zeta; -z, z^2], 'grid', {linspace(0, 1, 21), linspace(0, 1, 41)},...
 	'gridName', {'z', 'zeta'}, 'name', 'q');
@@ -168,8 +180,9 @@ tc.verifyEqual(qDiscrete2.imag.on(), -qDiscrete2Ctransp.imag.on());
 end % testCtranspose
 
 function testTranspose(tc)
+z = tc.TestData.sym.z;
+zeta = tc.TestData.sym.zeta;
 
-syms z zeta
 qSymbolic = quantity.Symbolic(...
 	[1+z*zeta, -zeta; -z, z^2], 'grid', {linspace(0, 1, 21), linspace(0, 1, 41)},...
 	'gridName', {'z', 'zeta'}, 'name', 'q');
@@ -183,7 +196,9 @@ tc.verifyEqual(qDiscrete(2,1).on(), qDiscreteTransp(1,2).on());
 end % testTranspose
 
 function testFlipGrid(tc)
-syms z zeta
+z = tc.TestData.sym.z;
+zeta = tc.TestData.sym.zeta;
+
 myGrid = linspace(0, 1, 11);
 f = quantity.Discrete(quantity.Symbolic([1+z+zeta; 2*zeta+sin(z)] + zeros(2, 1)*z*zeta, ...
 	'gridName', {'z', 'zeta'}, 'grid', {myGrid, myGrid}));
@@ -398,7 +413,8 @@ end
 
 function testDiag2Vec(testCase)
 % quantity.Symbolic
-syms z
+z = testCase.TestData.sym.z;
+
 do = quantity.Domain('z', linspace(0,1,3));
 myMatrixSymbolic = quantity.Symbolic([sin(0.5*z*pi)+1, 0; 0, 0.9-z/2], 'domain', do);
 myVectorSymbolic = diag2vec(myMatrixSymbolic);
@@ -463,7 +479,9 @@ testCase.verifyEqual(solution.on(), referenceResult1 + referenceResult2, 'AbsTol
 end
 
 function testSolveDVariableEqualQuantityNegative(testCase)
-syms z zeta
+z = testCase.TestData.sym.z;
+zeta = testCase.TestData.sym.zeta;
+
 assume(z>0 & z<1); assume(zeta>0 & zeta<1);
 myParameterGrid = linspace(0, 0.1, 5);
 Lambda = quantity.Symbolic(-0.1-z^2, ...%, -1.2+z^2]),...1+z*sin(z)
@@ -482,7 +500,8 @@ end
 
 function testSolveDVariableEqualQuantityComparedToSym(testCase)
 %% compare with symbolic implementation
-syms z
+z = testCase.TestData.sym.z;
+
 assume(z>0 & z<1);
 quanBSym = quantity.Symbolic(1+z, 'grid', {linspace(0, 1, 5)}, ...
 	'gridName', 'z', 'name', 'bSym', 'gridName', {'z'});
@@ -496,7 +515,8 @@ end
 
 function testSolveDVariableEqualQuantityAbsolut(testCase)
 %% compare with symbolic implementation
-syms z
+z = testCase.TestData.sym.z;
+
 assume(z>0 & z<1);
 quanBSym = quantity.Symbolic(1+z, 'grid', {linspace(0, 1, 11)}, ...
 	'gridName', 'z', 'name', 'bSym', 'gridName', {'z'});
@@ -505,7 +525,7 @@ quanBDiscrete = quantity.Discrete(quanBSym.on(), 'grid', {linspace(0, 1, 11)}, .
 
 solutionBDiscrete = quanBDiscrete.solveDVariableEqualQuantity();
 myGrid = solutionBDiscrete.grid{1};
-solutionBDiscreteDiff = solutionBDiscrete.diff(1, 's');
+solutionBDiscreteDiff = solutionBDiscrete.diff('s');
 quanOfSolutionOfS = zeros(length(myGrid), length(myGrid), 1);
 for icIdx = 1 : length(myGrid)
 	quanOfSolutionOfS(:, icIdx, :) = quanBSym.on(solutionBDiscrete.on({myGrid, myGrid(icIdx)}));
@@ -523,7 +543,8 @@ b = quantity.Discrete(ones(5, 1), 'grid', linspace(0, 1, 5), 'gridName', 'z', ..
 	'size', [1, 1], 'name', 'b');
 ab  = a*b;
 %
-syms z
+z = testCase.TestData.sym.z;
+
 c = quantity.Symbolic(1, 'grid', linspace(0, 1, 21), 'gridName', 'z', 'name', 'c');
 ac = a*c;
 
@@ -547,10 +568,8 @@ sinfun = quantity.Discrete(sin(z), 'grid', z, 'gridName', 'z');
 % very bad a the boundarys of the domain.
 Z = linspace(0.1, pi-0.1)';
 testCase.verifyTrue(numeric.near(sinfun.diff().on(Z), cos(Z), 1e-3));
-testCase.verifyTrue(numeric.near(sinfun.diff(2).on(Z), -sin(Z), 1e-3));
-testCase.verifyTrue(numeric.near(sinfun.diff(3).on(Z), -cos(Z), 1e-3));
-
-
+testCase.verifyTrue(numeric.near(sinfun.diff('z', 2).on(Z), -sin(Z), 1e-3));
+testCase.verifyTrue(numeric.near(sinfun.diff('z', 3).on(Z), -cos(Z), 1e-3));
 
 end
 
@@ -560,27 +579,26 @@ function testDiffConstant2d(testCase)
 myQuantity = quantity.Discrete(cat(3, 2*ones(11, 21), zNdgrid, zetaNdgrid), ...
 	'grid', {linspace(0, 1, 11), linspace(0, 1, 21)}, ...
 	'gridName', {'z', 'zeta'}, 'name', 'constant', 'size', [3, 1]);
-myQuantityDz = diff(myQuantity, 1, 'z');
-myQuantityDzeta = diff(myQuantity, 1, 'zeta');
-myQuantityDZzeta = diff(myQuantity, 1);
-myQuantityDZzeta2 = diff(myQuantity, 1, {'z', 'zeta'});
+myQuantityDz = diff(myQuantity, 'z', 1);
+myQuantityDzeta = diff(myQuantity, 'zeta', 1);
 
-testCase.verifyEqual(myQuantityDZzeta.on(), myQuantityDZzeta2.on());
+testCase.verifyError( @() diff(myQuantity), 'MATLAB:functionValidation:NotScalar' )
+testCase.verifyError( @() diff(myQuantity, 1, {'z', 'zeta'}), 'MATLAB:functionValidation:ClassConversionError' )
+testCase.verifyError( @() diff(myQuantity,  {'z', 'zeta'}), 'MATLAB:functionValidation:NotScalar' )
 
 % constant
 testCase.verifyEqual(myQuantityDz(1).on(), zeros(11, 21));
 testCase.verifyEqual(myQuantityDzeta(1).on(), zeros(11, 21));
-testCase.verifyEqual(myQuantityDZzeta(1).on(), zeros(11, 21));
+
 
 % zNdgrid
 testCase.verifyEqual(myQuantityDz(2).on(), ones(11, 21), 'AbsTol', 10*eps);
 testCase.verifyEqual(myQuantityDzeta(2).on(), zeros(11, 21), 'AbsTol', 10*eps);
-testCase.verifyEqual(myQuantityDZzeta(2).on(), zeros(11, 21), 'AbsTol', 10*eps);
 
 % zetaNdgrid
 testCase.verifyEqual(myQuantityDz(3).on(), zeros(11, 21), 'AbsTol', 10*eps);
 testCase.verifyEqual(myQuantityDzeta(3).on(), ones(11, 21), 'AbsTol', 10*eps);
-testCase.verifyEqual(myQuantityDZzeta(3).on(), zeros(11, 21), 'AbsTol', 10*eps);
+
 end
 
 function testOn(testCase)
@@ -669,7 +687,7 @@ zeta11 = quantity.Domain('zeta', z11_);
 a = quantity.Discrete(ones([z11.n, z11.n, 2, 2]), 'size', [2, 2], ...
 	'domain', [z11, zeta11], 'name', 'A');
 
-syms zeta
+zeta = testCase.TestData.sym.zeta;
 z10_ = linspace(0, 1, 10);
 
 b = quantity.Symbolic((eye(2, 2)), 'gridName', 'zeta', ...
@@ -683,7 +701,9 @@ end
 
 function testMTimesPointWise(tc)
 
-syms z zeta
+z = tc.TestData.sym.z;
+zeta = tc.TestData.sym.zeta;
+
 Z = linspace(0, pi, 51)';
 ZETA = linspace(0, pi / 2, 71)';
 
@@ -785,7 +805,8 @@ tGrid = linspace(pi, 1.1*pi, 51)';
 s = tGrid;
 [T, S] = ndgrid(tGrid, tGrid);
 
-syms sy t
+t = testCase.TestData.sym.t;
+sy = testCase.TestData.sym.sy;
 
 a = [ 1, sy; t, 1];
 b = [ sy; 2*sy];
@@ -842,7 +863,10 @@ testCase.verifyEqual( f.on(f(1).grid, f(1).gridName), V.on(f(1).grid, f(1).gridN
 tGrid = linspace(pi, 1.1*pi, 51)';
 s = tGrid;
 [T, S] = ndgrid(tGrid, tGrid);
-syms sy tt
+
+sy = testCase.TestData.sym.sy;
+t = testCase.TestData.sym.t;
+
 a = [ 1, sy; t, 1];
 
 A = zeros([size(T), size(a)]);
@@ -1287,12 +1311,12 @@ testCase.verifyEqual(quan10, shiftdim(quan.on({1, 0})));
 testCase.verifyEqual(quan01, shiftdim(quan.on({0, 1})));
 testCase.verifyEqual(quanZetaZ.on(), quan.on());
 testCase.verifyEqual(quant1.on(), squeeze(quan.on({quan(1).grid{1}, 1})));
-testCase.verifyEqual(quanZetaZ(1).gridName, {'zeta', 'z'})
-testCase.verifyEqual(quant1(1).gridName, {'t'})
+testCase.verifyEqual(quanZetaZ(1).gridName, ["zeta", "z"])
+testCase.verifyEqual(quant1(1).gridName, ["t"])
 
 quanZetaZetaReference = misc.diagNd(quan.on({linspace(0, 1, 41), linspace(0, 1, 41)}), [1, 2]);
 testCase.verifyEqual(quanZetaZeta.on(), quanZetaZetaReference);
-testCase.verifyEqual(quanEta(1).gridName, {'eta', 'zeta'})
+testCase.verifyEqual(quanEta(1).gridName, ["eta", "zeta"])
 testCase.verifyEqual(quanPt.on(), shiftdim(quan.on({1, quan(1).grid{2}})));
 
 myQuantity = quanZetaZeta.subs(quantity.Domain('zeta', linspace(0,1,3)));
diff --git a/+unittests/+quantity/testDomain.m b/+unittests/+quantity/testDomain.m
index 7fff80101e082291c5f6dc6184b690daf78f08d5..444efef2aaa5fad3b943058f75ee941cd66c567e 100644
--- a/+unittests/+quantity/testDomain.m
+++ b/+unittests/+quantity/testDomain.m
@@ -29,9 +29,9 @@ d = testCase.TestData.d;
 testCase.verifyEqual( d.find('a'), a)
 testCase.verifyEqual( d.find('b', 'c'), [b c])
 testCase.verifyEqual( d.find({'b', 'c'}), [b c])
-testCase.verifyEqual( d.find(b), b);
-testCase.verifyEqual( d.find({a c}), [a c])
-testCase.verifyEqual( d.find([a b]), [a b])
+testCase.verifyEqual( d.find(b.name), b);
+testCase.verifyEqual( d.find({a.name c.name}), [a c])
+testCase.verifyEqual( d.find([a.name b.name]), [a b])
 
 end
 
@@ -45,7 +45,7 @@ c = quantity.Domain('c', z);
 d = [a b c];
 
 testCase.verifyEqual(d.find('a'), a);
-testCase.verifyEqual(d.find(b), b);
+testCase.verifyEqual(d.find(b.name), b);
 testCase.verifyEqual(d.find('c'), c);
 
 end
@@ -88,7 +88,7 @@ d = [a b c];
 [d__, I__] = d.sort('ascend');
 d___ = d.sort();
 
-testCase.verifyEqual({d_.name}, {'c', 'b', 'a'});
+testCase.verifyEqual([d_.name], ["c", "b", "a"]);
 testCase.verifyEqual({d__.name}, {d.name});
 testCase.verifyEqual({d__.name}, {d___.name});
 testCase.verifyEqual( I_, flip(I__) );
@@ -168,16 +168,22 @@ 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.name], ["p", "s", "t", "z"]);
 testCase.verifyEqual({joinedDomainAB.grid}, {s(:), s(:), t(:), z(:)});
-testCase.verifyEqual({joinedDomainCC.name}, {'p'});
+testCase.verifyEqual(joinedDomainCC.name, "p");
 testCase.verifyEqual({joinedDomainCC.grid}, {t(:)});
-testCase.verifyEqual({joinedDomainBC.name}, {'p', 's'});
+testCase.verifyEqual([joinedDomainBC.name], ["p", "s"]);
 testCase.verifyEqual({joinedDomainBC.grid}, {s(:), t(:)});
 testCase.verifyEqual({joinedDomain.grid}, {t(:), z(:)});
+
+d1 = [Z Z T];
+d2 = quantity.Domain.empty;
+
+testCase.verifyEqual( d1.join( d2 ), [Z T] )
+
 end
 
-function testGridIndex(testCase)
+function testFind_Index(testCase)
 z = linspace(0, 2*pi, 71)';
 t = linspace(0, 3*pi, 51);
 s = linspace(0, 1);
diff --git a/+unittests/+quantity/testSymbolic.m b/+unittests/+quantity/testSymbolic.m
index 96313753a8d395655240d3916b8c7b0d3c77f846..8a0e2ae64b002308bb57b02c77676e831ff81a2b 100644
--- a/+unittests/+quantity/testSymbolic.m
+++ b/+unittests/+quantity/testSymbolic.m
@@ -140,13 +140,8 @@ fdz = quantity.Symbolic([zeta, 1; 0, 0], ...
 fdzeta = quantity.Symbolic([z, 0; 1, 0], ...
 	'gridName', {'z', 'zeta'}, 'grid', {myGrid, myGrid});
 
-fRefTotal = 0*f.on();
-fRefTotal(:,:,1,1) = 1;
-testCase.verifyEqual(on(diff(f)), fRefTotal);
-testCase.verifyEqual(on(diff(f, 1, {'z', 'zeta'})), fRefTotal);
-testCase.verifyEqual(on(diff(f, 2, {'z', 'zeta'})), 0*fRefTotal);
-testCase.verifyEqual(on(diff(f, 1, 'z')), fdz.on());
-testCase.verifyEqual(on(diff(f, 1, 'zeta')), fdzeta.on());
+testCase.verifyEqual(on(diff(f, 'z')), fdz.on());
+testCase.verifyEqual(on(diff(f, 'zeta')), fdzeta.on());
 
 end
 
@@ -483,11 +478,10 @@ f2 = cosh(z * pi);
 f = quantity.Symbolic([f1 ; f2 ], 'name', 'f', 'domain', ...
 	quantity.Domain('z', linspace(0,1,11)));
 
-d2f = f.diff(2);
-d1f = f.diff(1);
+d2f = f.diff("z", 2);
+d1f = f.diff("z", 1);
 
 %%
-
 F1 = symfun([f1; f2], z);
 D1F = symfun([diff(f1,1); diff(f2,1)], z);
 D2F = symfun([diff(f1,2); diff(f2,2)], z);
@@ -500,8 +494,8 @@ R3 = D2F(z);
 
 %%
 verifyTrue(testCase, numeric.near(double([R1{:}]), f.on(), 1e-12));
-verifyTrue(testCase, numeric.near(double([R2{:}]), f.diff(1).on(), 1e-12));
-verifyTrue(testCase, numeric.near(double([R3{:}]), f.diff(2).on(), 1e-12));
+verifyTrue(testCase, numeric.near(double([R2{:}]), f.diff('z', 1).on(), 1e-12));
+verifyTrue(testCase, numeric.near(double([R3{:}]), f.diff('z', 2).on(), 1e-12));
 
 end
 
@@ -683,7 +677,7 @@ Q_var = quantity.Symbolic(z, 'domain', quantity.Domain('z', 1:2));
 
 testCase.verifyError(...
 	@() quantity.Symbolic(z, 'domain', quantity.Domain.empty()), ...
-	'quantity:Symbolic:Initialisation' );
+	'symbolic:sym:matlabFunction:InvalidVars' );
 
 testCase.verifyTrue(Q_constant.isConstant);
 testCase.verifyFalse(Q_var.isConstant);
@@ -809,16 +803,16 @@ assume(z>0 & z<1); assume(zeta>0 & zeta<1);
 F = quantity.Symbolic(zeros(1), 'grid', ...
 	{linspace(0, 1, 3), linspace(0, 1, 3)}, 'gridName', {'z', 'zeta'});
 Feta = F.subs('z', 'eta');
-testCase.verifyEqual(Feta(1).gridName, {'eta', 'zeta'});
+testCase.verifyEqual(Feta(1).gridName, ["eta", "zeta"]);
 
 %%
 fMessy = quantity.Symbolic(x^1*y^2*z^4*w^5, 'domain', [X, Y, Z, W]);
 fMessyA = fMessy.subs({'x', 'y', 'z', 'w'}, {'a', 'a', 'a', 'a'});
 fMessyYY = fMessy.subs({'w', 'x', 'z'}, {'xi', 'y', 'y'});
 
-testCase.verifyEqual(fMessyA.gridName, {'a'});
+testCase.verifyEqual(fMessyA.gridName, "a");
 testCase.verifyEqual(fMessyA.grid, {linspace(0, 1, 11)'});
-testCase.verifyEqual(fMessyYY.gridName, {'y', 'xi'});
+testCase.verifyEqual(fMessyYY.gridName, ["y", "xi"]);
 testCase.verifyEqual(fMessyYY.grid, {linspace(0, 1, 11)', linspace(0, 1, 3)'});
 % % sub multiple numerics -> not implemented yet
 % f11 = f.subs('x', [1; 2]);