--- a/+blockmatrix/toMatrix.m	Fri Oct 20 16:01:43 2017 +0200
+++ b/+blockmatrix/toMatrix.m	Fri Oct 20 23:31:20 2017 +0200
@@ -12,16 +12,12 @@
 
     A = sparse(N,M);
 
-    n_ind = [0 cumsum(n)];
-    m_ind = [0 cumsum(m)];
-
     for i = 1:size(bm,1)
         for j = 1:size(bm,2)
             if isempty(bm{i,j})
-                continue
+                bm{i,j} = sparse(n(i),m(j));
             end
-            % TODO: If this ever fails for large matrices. Try cell2mat instead.
-            A(n_ind(i)+1:n_ind(i+1),m_ind(j)+1:m_ind(j+1)) = bm{i,j};
         end
     end
+    A = cell2mat(bm);
 end

--- a/+grid/evalOn.m	Fri Oct 20 16:01:43 2017 +0200
+++ b/+grid/evalOn.m	Fri Oct 20 23:31:20 2017 +0200
@@ -7,42 +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);
 
-
-    % Get coordinates and convert to cell array for easier use as a parameter
-    x = num2cell(g.points());
+    gf = func(x{:});
+    gf = reorderComponents(gf, k);
+end
 
-    % Find the number of components
-    if size(x,1) ~= 0
-        x0 = x(1,:);
-    else
-        x0 = num2cell(ones(1,size(x,2)));
-    end
+% 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
 
-    if size(f0,2) ~= 1
-        error('grid:evalOn:VectorValuedWrongDim', 'A vector valued function must be given as a column vector')
+% 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
-
-    gf = zeros(g.N*k, 1);
-    for i = 1:g.N
-        % (1 + (i-1)*k):(i*k)
-        % func(x{i,:})
-        gf((1 + (i-1)*k):(i*k)) = func(x{i,:});
-    end
-end
\ No newline at end of file
+end

--- a/+grid/evalOnTest.m	Fri Oct 20 16:01:43 2017 +0200
+++ b/+grid/evalOnTest.m	Fri Oct 20 23:31:20 2017 +0200
@@ -31,7 +31,7 @@
     cases = {
         {getTestGrid('1d'), @(x,y)x-y},
         {getTestGrid('2d'), @(x)x    },
-    }
+    };
 
     for i = 1:length(cases)
         g = cases{i}{1};
@@ -111,9 +111,9 @@
 
 
 function testInputErrorVectorValued(testCase)
-     in  = {
+    in  = {
         [1,2,3],
-        @(x,y)[x,-y];
+        @(x,y)[x,-y],
     };
 
     g = getTestGrid('2d');