--- a/+noname/calculateErrors.m	Wed Jul 25 18:53:07 2018 -0700
+++ b/+noname/calculateErrors.m	Thu Jul 26 10:53:05 2018 -0700
@@ -4,27 +4,40 @@
 % m are grid size parameters.
 % N are number of timesteps to use for each gird size
 % timeOpt are options for the timeStepper
+% errorFun is a function_handle taking 2 or 3 arguments, errorFun(trueSolution, approxSolution), errorFun(trueSolution, approxSolution, discr)
 function e = calculateErrors(schemeFactory, T, m, N, errorFun, timeOpt)
+    %TODO: Ability to choose paralell or not
     assertType(schemeFactory, 'function_handle');
     assertNumberOfArguments(schemeFactory, 1);
     assertScalar(T);
     assert(length(m) == length(N), 'Vectors m and N must have the same length');
     assertType(errorFun, 'function_handle');
-    assertNumberOfArguments(errorFun, 2);
-    default_arg('timeOpt');
+
+    if ~ismember(nargin(errorFun), [2,3])
+        error('sbplib:noname:calculateErrors:wrongNumberOfArguments', '"%s" must have 2 or 3, found %d', toString(errorFun), nargin(errorFun));
+    end
 
-    e = [];
-    for i = 1:length(m)
+    default_arg('timeOpt', struct());
+
+
+    e = zeros(1,length(m));
+    parfor i = 1:length(m)
         done = timeTask('m = %3d ', m(i));
 
         [discr, trueSolution] = schemeFactory(m(i));
 
-        timeOpt.k = T/N(i);
-        ts = discr.getTimestepper(timeOpt);
+        timeOptTemp = timeOpt;
+        timeOptTemp.k = T/N(i);
+        ts = discr.getTimestepper(timeOptTemp);
         ts.stepTo(N(i), true);
         approxSolution = discr.getTimeSnapshot(ts);
 
-        e(i) = errorFun(trueSolution, approxSolution);
+        switch nargin(errorFun)
+            case 2
+                e(i) = errorFun(trueSolution, approxSolution);
+            case 3
+                e(i) = errorFun(trueSolution, approxSolution, discr);
+        end
 
         fprintf('e = %.4e', e(i))
         done()