--- a/+grid/evalOn.m	Tue May 22 13:29:47 2018 -0700
+++ b/+grid/evalOn.m	Tue Jul 24 20:14:29 2018 -0700
@@ -7,68 +7,34 @@
 function gf = evalOn(g, func)
     if ~isa(func, 'function_handle')
         % We should have a constant.
-        if size(func,2) ~= 1
-            error('grid:evalOn:VectorValuedWrongDim', 'A vector valued function must be given as a column vector')
-        end
+        assert(size(func,2) == 1,'grid:evalOn:VectorValuedWrongDim', 'A vector valued function must be given as a column vector');
 
         gf = repmat(func,[g.N, 1]);
         return
     end
     % func should now be a function_handle
+    assert(g.D == nargin(func),'grid:evalOn:WrongNumberOfInputs', 'The number of inputs of the function must match the dimension of the domain.')
 
-    if g.D ~= nargin(func)
-        g.D
-        nargin(func)
-        error('grid:evalOn:WrongNumberOfInputs', 'The number of inputs of the function must match the dimension of the domain.')
-    end
+    x = num2cell(g.points(),1);
+    k = numberOfComponents(func);
 
+    gf = func(x{:});
+    gf = reorderComponents(gf, k);
+end
 
-    % Get coordinates
-    x = g.points();
+% Find the number of vector components of func
+function k = numberOfComponents(func)
+    x0 = num2cell(ones(1,nargin(func)));
+    f0 = func(x0{:});
+    assert(size(f0,2) == 1, 'grid:evalOn:VectorValuedWrongDim', 'A vector valued function must be given as a column vector');
+    k = length(f0);
+end
 
-    % Find the number of components
-    if size(x,1) ~= 0
-        x0 = x(1,:);
-    else
-        x0 = num2cell(ones(1,size(x,2)));
+% Reorder the components of the function to sit together
+function gf = reorderComponents(a, k)
+    N = length(a)/k;
+    gf = zeros(N*k, 1);
+    for i = 1:k
+        gf(i:k:end) = a((i-1)*N + 1 : i*N);
     end
-    
-    dim = length(x0);
-    % Evaluate f0 = func(x0(1),x0(2),...,x0(dim));
-    if(dim == 1)
-        f0 = func(x0);
-    else
-        eval_str = 'f0 = func(x0(1)';
-        for i = 2:dim
-            eval_str = [eval_str, sprintf(',x0(%d)',i)];
-        end
-        eval_str = [eval_str, ');'];
-        eval(eval_str);
-    end
-
-    % k = number of components
-    k = length(f0);
-
-    if size(f0,2) ~= 1
-        error('grid:evalOn:VectorValuedWrongDim', 'A vector valued function must be given as a column vector')
-    end
-
-    % Evaluate gf = func(x(:,1),x(:,2),...,x(:,dim));
-    if(dim == 1)
-        gf = func(x);
-    else
-        eval_str = 'gf = func(x(:,1)';
-        for i = 2:dim
-            eval_str = [eval_str, sprintf(',x(:,%d)',i)];
-        end
-        eval_str = [eval_str, ');'];
-        eval(eval_str);
-    end
-    
-    % Reorganize gf
-    gf_temp = gf;
-    gf = zeros(g.N*k, 1);
-    for i = 1:k
-        gf(i:k:end) = gf_temp((i-1)*g.N + 1 : i*g.N);
-    end
-end
\ No newline at end of file
+end