--- a/+anim/setup_time_quantity_plot.m	Tue Jul 24 19:59:54 2018 -0700
+++ b/+anim/setup_time_quantity_plot.m	Tue Jul 24 20:01:29 2018 -0700
@@ -16,7 +16,7 @@
         if ishandle(axis_handle)
             % t = [t t_now];
             for j = 1:length(yfun)
-                addpoints(plot_handles(j),t_now,yfun{j}(varargin{:}));
+                addpoints(plot_handles(j),t_now,full(yfun{j}(varargin{:})));
             end
 
             [t,~] = getpoints(plot_handles(1));

--- a/+blockmatrix/toMatrix.m	Tue Jul 24 19:59:54 2018 -0700
+++ b/+blockmatrix/toMatrix.m	Tue Jul 24 20:01:29 2018 -0700
@@ -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/+draw/prompt_point.m	Tue Jul 24 19:59:54 2018 -0700
+++ b/+draw/prompt_point.m	Tue Jul 24 20:01:29 2018 -0700
@@ -1,22 +1,23 @@
-function [p, button] = prompt_point(s,varargin)
+function [p, button] = prompt_point(s, varargin)
     default_arg('s',[])
 
     set(gcf,'Pointer','crosshair')
 
     if ~isempty(s)
-        fprintf(s,varargin{:});
+        fprintf(s, varargin{:});
     end
 
-    a = gca;
+    fh = gcf();
+    ah = gca();
 
-    function get_point(src,event)
-        cp = a.CurrentPoint;
+    function get_point(src, event)
+        cp = ah.CurrentPoint;
         p = cp(1,1:2)';
-        a.ButtonDownFcn = [];
+        fh.WindowButtonUpFcn = [];
     end
 
-    a.ButtonDownFcn = @get_point;
-    waitfor(a,'ButtonDownFcn', [])
+    fh.WindowButtonUpFcn = @get_point;
+    waitfor(fh,'WindowButtonUpFcn', [])
 
     set(gcf,'Pointer','arrow')
 

--- a/+grid/Cartesian.m	Tue Jul 24 19:59:54 2018 -0700
+++ b/+grid/Cartesian.m	Tue Jul 24 20:01:29 2018 -0700
@@ -15,7 +15,7 @@
             obj.d = length(varargin);
 
             for i = 1:obj.d
-                assert(isvector(varargin{i}), 'Coordinate inputs must be a vectors.')
+                assert(isnumeric(varargin{i}), 'Coordinate inputs must be vectors.')
 
                 obj.x{i} = varargin{i};
                 obj.m(i) = length(varargin{i});

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+grid/TODO.txt	Tue Jul 24 20:01:29 2018 -0700
@@ -0,0 +1,1 @@
+% TODO: Rename grid package. name conflicts with built in function

--- a/+grid/evalOn.m	Tue Jul 24 19:59:54 2018 -0700
+++ b/+grid/evalOn.m	Tue Jul 24 20:01:29 2018 -0700
@@ -7,43 +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);
-    % keyboard
-    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	Tue Jul 24 19:59:54 2018 -0700
+++ b/+grid/evalOnTest.m	Tue Jul 24 20:01:29 2018 -0700
@@ -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');

--- a/+multiblock/DefCurvilinear.m	Tue Jul 24 19:59:54 2018 -0700
+++ b/+multiblock/DefCurvilinear.m	Tue Jul 24 20:01:29 2018 -0700
@@ -48,13 +48,14 @@
             g = multiblock.Grid(grids, obj.connections, obj.boundaryGroups);
         end
 
-        function show(obj, label, gridLines, varargin)
+        function h = show(obj, label, gridLines, varargin)
             default_arg('label', 'name')
             default_arg('gridLines', false);
 
+            h = [];
             if isempty('label') && ~gridLines
                 for i = 1:obj.nBlocks
-                    obj.blockMaps{i}.show(2,2);
+                    h = [h, obj.blockMaps{i}.show(2,2)];
                 end
                 axis equal
                 return
@@ -63,7 +64,7 @@
             if gridLines
                 ms = obj.getGridSizes(varargin{:});
                 for i = 1:obj.nBlocks
-                    obj.blockMaps{i}.show(ms{i}(1),ms{i}(2));
+                    h = [h, obj.blockMaps{i}.show(ms{i}(1),ms{i}(2))];
                 end
             end
 
@@ -76,7 +77,7 @@
                     for i = 1:obj.nBlocks
                         labels{i} = num2str(i);
                     end
-                otherwise
+                case 'none'
                     axis equal
                     return
             end

--- a/+multiblock/DiffOp.m	Tue Jul 24 19:59:54 2018 -0700
+++ b/+multiblock/DiffOp.m	Tue Jul 24 20:01:29 2018 -0700
@@ -10,13 +10,13 @@
     end
 
     methods
-        function obj = DiffOp(doHand, grid, order, doParam)
+        function obj = DiffOp(doHand, g, order, doParam)
             %  doHand -- may either be a function handle or a cell array of
             %            function handles for each grid. The function handle(s)
             %            should be on the form do = doHand(grid, order, ...)
             %            Additional parameters for each doHand may be provided in
             %            the doParam input.
-            %    grid -- a multiblock grid
+            %       g -- a multiblock grid
             %   order -- integer specifying the order of accuracy
             % doParam -- may either be a cell array or a cell array of cell arrays
             %            for each block. If it is a cell array with length equal
@@ -26,9 +26,9 @@
             %            extra parameters to all doHand: doHand(..., doParam{:})
             default_arg('doParam', [])
 
-            [getHand, getParam] = parseInput(doHand, grid, doParam);
+            [getHand, getParam] = parseInput(doHand, g, doParam);
 
-            nBlocks = grid.nBlocks();
+            nBlocks = g.nBlocks();
 
             obj.order = order;
 
@@ -40,7 +40,7 @@
                 if ~iscell(p)
                     p = {p};
                 end
-                obj.diffOps{i} = h(grid.grids{i}, order, p{:});
+                obj.diffOps{i} = h(g.grids{i}, order, p{:});
             end
 
 
@@ -53,7 +53,7 @@
 
 
             % Build the differentiation matrix
-            obj.blockmatrixDiv = {grid.Ns, grid.Ns};
+            obj.blockmatrixDiv = {g.Ns, g.Ns};
             D = blockmatrix.zero(obj.blockmatrixDiv);
             for i = 1:nBlocks
                 D{i,i} = obj.diffOps{i}.D;
@@ -61,7 +61,7 @@
 
             for i = 1:nBlocks
                 for j = 1:nBlocks
-                    intf = grid.connections{i,j};
+                    intf = g.connections{i,j};
                     if isempty(intf)
                         continue
                     end
@@ -77,14 +77,15 @@
                 end
             end
             obj.D = blockmatrix.toMatrix(D);
+            obj.grid = g;
 
 
-            function [getHand, getParam] = parseInput(doHand, grid, doParam)
-                if ~isa(grid, 'multiblock.Grid')
+            function [getHand, getParam] = parseInput(doHand, g, doParam)
+                if ~isa(g, 'multiblock.Grid')
                     error('multiblock:DiffOp:DiffOp:InvalidGrid', 'Requires a multiblock grid.');
                 end
 
-                if iscell(doHand) && length(doHand) == grid.nBlocks()
+                if iscell(doHand) && length(doHand) == g.nBlocks()
                     getHand = @(i)doHand{i};
                 elseif isa(doHand, 'function_handle')
                     getHand = @(i)doHand;
@@ -104,7 +105,7 @@
 
                 % doParam is a non-empty cell-array
 
-                if length(doParam) == grid.nBlocks() && all(cellfun(@iscell, doParam))
+                if length(doParam) == g.nBlocks() && all(cellfun(@iscell, doParam))
                     % doParam is a cell-array of cell-arrays
                     getParam = @(i)doParam{i};
                     return
@@ -116,7 +117,7 @@
 
         function ops = splitOp(obj, op)
             % Splits a matrix operator into a cell-matrix of matrix operators for
-            % each grid.
+            % each g.
             ops = sparse2cell(op, obj.NNN);
         end
 
@@ -143,6 +144,27 @@
             end
         end
 
+        function op = getBoundaryQuadrature(obj, boundary)
+            opName = 'H';
+            switch class(boundary)
+                case 'cell'
+                    localOpName = [opName '_' boundary{2}];
+                    blockId = boundary{1};
+                    op = obj.diffOps{blockId}.(localOpName);
+
+                    return
+                case 'multiblock.BoundaryGroup'
+                    N = length(boundary);
+                    H_bm = cell(N,N);
+                    for i = 1:N
+                        H_bm{i,i} = obj.getBoundaryQuadrature(boundary{i});
+                    end
+                    op = blockmatrix.toMatrix(H_bm);
+                otherwise
+                    error('Unknown boundary indentifier')
+            end
+        end
+
         % Creates the closure and penalty matrix for a given boundary condition,
         %    boundary -- the name of the boundary on the form {id,name} where
         %                id is the number of a block and name is the name of a
@@ -194,6 +216,7 @@
                 p{I} = blockPenalty;
                 penalty = blockmatrix.toMatrix(p);
             else
+                % TODO: used by beam equation, should be eliminated. SHould only set one BC per call
                 for i = 1:length(blockPenalty)
                     div{2} = size(blockPenalty{i}, 2); % Penalty is a column vector
                     p = blockmatrix.zero(div);

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+multiblock/Laplace.m	Tue Jul 24 20:01:29 2018 -0700
@@ -0,0 +1,56 @@
+classdef Laplace < scheme.Scheme
+    properties
+        grid
+        order
+        mbDiffOp
+
+        D
+        H
+        J
+    end
+    methods
+        function obj = Laplace(g, order, a, b, opGen)
+            default_arg('order', 4);
+            default_arg('a', 1);
+            default_arg('b', 1);
+            default_arg('opGen', @sbp.D4Variable);
+
+            obj.grid = g;
+            obj.order = order;
+            obj.mbDiffOp = multiblock.DiffOp(@scheme.LaplaceCurvilinear, obj.grid, order, {a,b,opGen});
+
+            obj.D = obj.mbDiffOp.D;
+            obj.J = obj.jacobian();
+            obj.H = obj.mbDiffOp.H;
+        end
+
+        function s = size(obj)
+            s = size(obj.mbDiffOp);
+        end
+
+        function J = jacobian(obj)
+            N = obj.grid.nBlocks;
+            J = cell(N,N);
+
+            for i = 1:N
+                J{i,i} = obj.mbDiffOp.diffOps{i}.J;
+            end
+            J = blockmatrix.toMatrix(J);
+        end
+
+        function op = getBoundaryOperator(obj, opName, boundary)
+            op = getBoundaryOperator(obj.mbDiffOp, opName, boundary);
+        end
+
+        function op = getBoundaryQuadrature(obj, boundary)
+            op = getBoundaryQuadrature(obj.mbDiffOp, boundary);
+        end
+
+        function [closure, penalty] = boundary_condition(obj,boundary,type) % TODO: Change name to boundaryCondition
+            [closure, penalty] = boundary_condition(obj.mbDiffOp, boundary, type);
+        end
+        function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
+            error('Not implemented')
+        end
+    end
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+noname/Animation.m	Tue Jul 24 20:01:29 2018 -0700
@@ -0,0 +1,75 @@
+classdef Animation < handle
+    properties
+        timeStepper
+        representationMaker
+        updaters
+    end
+
+    % add input validation
+
+    methods
+        function obj = Animation(timeStepper, representationMaker, updaters);
+            obj.timeStepper = timeStepper;
+            obj.updaters = updaters;
+            obj.representationMaker = representationMaker;
+        end
+
+        function update(obj, r)
+            for i = 1:length(obj.updaters)
+                obj.updaters{i}(r);
+            end
+            drawnow
+        end
+
+        function run(obj, tEnd, timeModifier, do_pause)
+            default_arg('do_pause', false)
+
+            function next_t = G(next_t)
+                obj.timeStepper.evolve(next_t);
+                r = obj.representationMaker(obj.timeStepper);
+                obj.update(r);
+
+                if do_pause
+                    pause
+                end
+            end
+
+            anim.animate(@G, obj.timeStepper.t, tEnd, timeModifier);
+        end
+
+        function step(obj, tEnd, do_pause)
+            default_arg('do_pause', false)
+
+            while obj.timeStepper.t < tEnd
+                obj.timeStepper.step();
+
+                r = obj.representationMaker(obj.timeStepper);
+                obj.update(r);
+
+                % TODO: Make it never go faster than a certain fram rate
+
+                if do_pause
+                    pause
+                end
+            end
+        end
+
+        function saveMovie(obj, tEnd, timeModifier, figureHandle, dirname)
+            save_frame = anim.setup_fig_mov(figureHandle, dirname);
+
+            function next_t = G(next_t)
+                obj.timeStepper.evolve(next_t);
+                r = obj.representationMaker(obj.timeStepper);
+                obj.update(r);
+
+                save_frame();
+            end
+
+            fprintf('Generating and saving frames to: ..\n')
+            anim.animate(@G, obj.timeStepper.t, tEnd, timeModifier);
+            fprintf('Generating movies...\n')
+            cmd = sprintf('bash %s/+anim/make_movie.sh %s', sbplibLocation(),dirname);
+            system(cmd);
+        end
+    end
+end

--- a/+noname/animate.m	Tue Jul 24 19:59:54 2018 -0700
+++ b/+noname/animate.m	Tue Jul 24 20:01:29 2018 -0700
@@ -1,10 +1,10 @@
-% hand = noname.animate(discretization, time_modifier, Tend, dirname, opt)
+% noname.animate(discretization, time_modifier, Tend, dirname, opt)
 %
 % Example:
 %      noname.animate(discr,timemodifier,tend)
 %      noname.animate(discr,1, [tstart tend],'my_mov', opt)
 
-function hand = animate(discretization, time_modifier, Tend, dirname, opt)
+function animate(discretization, time_modifier, Tend, dirname, opt)
     default_arg('time_modifier', 1);
     default_arg('Tend', Inf);
     default_arg('dirname', '');
@@ -87,7 +87,8 @@
         pause
         anim.animate(@G, Tstart, Tend, time_modifier);
     else
-        while true
+        pause
+        while ts.t < Tend
             ts.step();
             sol = discretization.getTimeSnapshot(ts);
             update(sol);

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+noname/calculateErrors.m	Tue Jul 24 20:01:29 2018 -0700
@@ -0,0 +1,53 @@
+% [discr, trueSolution] =  schemeFactory(m)
+%     where trueSolution should be a timeSnapshot of the true solution a time T
+% T is the end time
+% 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');
+
+    if ~ismember(nargin(errorFun), [2,3])
+        error('sbplib:noname:calculateErrors:wrongNumberOfArguments', '"%s" must have 2 or 3, found %d', toString(errorFun), nargin(errorFun));
+    end
+
+    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));
+
+        timeOptTemp = timeOpt;
+        timeOptTemp.k = T/N(i);
+        ts = discr.getTimestepper(timeOptTemp);
+        ts.stepTo(N(i), true);
+        approxSolution = discr.getTimeSnapshot(ts);
+
+        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()
+    end
+    fprintf('\n')
+end
+
+
+%% Example error function
+% u_true = grid.evalOn(dr.grid, @(x,y)trueSolution(T,x,y));
+% err = u_true-u_false;
+% e(i) = norm(err)/norm(u_true);
+% % e(i) = sqrt(err'*d.H*d.J*err/(u_true'*d.H*d.J*u_true));

--- a/+parametrization/Curve.m	Tue Jul 24 19:59:54 2018 -0700
+++ b/+parametrization/Curve.m	Tue Jul 24 20:01:29 2018 -0700
@@ -181,8 +181,8 @@
         end
 
         function D = mirror(C, a, b)
-            assert_size(a,[2,1]);
-            assert_size(b,[2,1]);
+            assertSize(a,[2,1]);
+            assertSize(b,[2,1]);
 
             g = C.g;
             gp = C.gp;
@@ -219,8 +219,8 @@
         end
 
         function D = rotate(C,a,rad)
-            assert_size(a, [2,1]);
-            assert_size(rad, [1,1]);
+            assertSize(a, [2,1]);
+            assertSize(rad, [1,1]);
             g = C.g;
             gp = C.gp;
 

--- a/+parametrization/Ti.m	Tue Jul 24 19:59:54 2018 -0700
+++ b/+parametrization/Ti.m	Tue Jul 24 20:01:29 2018 -0700
@@ -21,16 +21,29 @@
             D = g4(0);
 
             function o = S_fun(u,v)
+                if isrow(u) && isrow(v)
+                    flipped = false;
+                else
+                    flipped = true;
+                    u = u';
+                    v = v';
+                end
+
                 x1 = g1(u);
                 x2 = g2(v);
                 x3 = g3(1-u);
                 x4 = g4(1-v);
+
                 o1 = (1-v).*x1(1,:) + u.*x2(1,:) + v.*x3(1,:) + (1-u).*x4(1,:) ...
-                    -((1-u)*(1-v).*A(1,:) + u*(1-v).*B(1,:) + u*v.*C(1,:) + (1-u)*v.*D(1,:));
+                    -((1-u).*(1-v).*A(1,:) + u.*(1-v).*B(1,:) + u.*v.*C(1,:) + (1-u).*v.*D(1,:));
                 o2 = (1-v).*x1(2,:) + u.*x2(2,:) + v.*x3(2,:) + (1-u).*x4(2,:) ...
-                    -((1-u)*(1-v).*A(2,:) + u*(1-v).*B(2,:) + u*v.*C(2,:) + (1-u)*v.*D(2,:));
+                    -((1-u).*(1-v).*A(2,:) + u.*(1-v).*B(2,:) + u.*v.*C(2,:) + (1-u).*v.*D(2,:));
 
-                o = [o1;o2];
+                if ~flipped
+                    o = [o1;o2];
+                else
+                    o = [o1'; o2'];
+                end
             end
 
             obj.S = @S_fun;
@@ -116,13 +129,13 @@
             S = obj.S;
 
             if(nu>2 || nv>2)
-                h_grid = obj.plot(nu,nv);
-                set(h_grid,'Color',[0 0.4470 0.7410]);
+                h.grid = obj.plot(nu,nv);
+                set(h.grid,'Color',[0 0.4470 0.7410]);
             end
 
-            h_bord = obj.plot(2,2);
-            set(h_bord,'Color',[0.8500 0.3250 0.0980]);
-            set(h_bord,'LineWidth',2);
+            h.border = obj.plot(2,2);
+            set(h.border,'Color',[0.8500 0.3250 0.0980]);
+            set(h.border,'LineWidth',2);
         end
 
 

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+parametrization/TiTest.m	Tue Jul 24 20:01:29 2018 -0700
@@ -0,0 +1,52 @@
+function tests = TiTest()
+    tests = functiontests(localfunctions);
+end
+
+function testScalarInput(testCase)
+    ti = getMinimumTi();
+
+    cases = {
+        % {u, v, out},
+        {0, 0, [1; 2]},
+        {0, 1, [1; 4]},
+        {1, 0, [3; 2]},
+        {1, 1, [3; 4]},
+        {0.5, 0.5, [2; 3]},
+    };
+
+    for i = 1:length(cases)
+        u = cases{i}{1};
+        v = cases{i}{2};
+        expected = cases{i}{3};
+
+        testCase.verifyEqual(ti.S(u,v), expected, sprintf('Case: %d',i));
+    end
+end
+
+function testRowVectorInput(testCase)
+    ti = getMinimumTi();
+
+    u = [0, 0.5, 1];
+    v = [0, 0, 0.5];
+    expected = [
+        1, 2, 3;
+        2, 2, 3;
+    ];
+
+    testCase.verifyEqual(ti.S(u,v), expected);
+end
+
+function testColumnvectorInput(testCase)
+   ti = getMinimumTi();
+
+    u = [0; 0.5; 1];
+    v = [0; 0; 0.5];
+    expected = [1; 2; 3; 2; 2; 3];
+
+    testCase.verifyEqual(ti.S(u,v), expected);
+end
+
+
+function ti = getMinimumTi()
+    ti = parametrization.Ti.rectangle([1; 2], [3; 4]);
+end
\ No newline at end of file

--- a/+scheme/Beam.m	Tue Jul 24 19:59:54 2018 -0700
+++ b/+scheme/Beam.m	Tue Jul 24 20:01:29 2018 -0700
@@ -126,6 +126,44 @@
                     penalty{1} = -obj.Hi*tau;
                     penalty{1} = -obj.Hi*sig;
 
+                case 'e'
+                    alpha = obj.alpha;
+                    tuning = 1.1;
+
+                    tau1 = tuning * alpha/delt;
+                    tau4 = s*alpha;
+
+                    tau = tau1*e+tau4*d3;
+
+                    closure = obj.Hi*tau*e';
+                    penalty = -obj.Hi*tau;
+                case 'd1'
+                    alpha = obj.alpha;
+
+                    tuning = 1.1;
+
+                    sig2 = tuning * alpha/gamm;
+                    sig3 = -s*alpha;
+
+                    sig = sig2*d1+sig3*d2;
+
+                    closure = obj.Hi*sig*d1';
+                    penalty = -obj.Hi*sig;
+
+                case 'd2'
+                    a = obj.alpha;
+
+                    tau =  s*a*d1;
+
+                    closure = obj.Hi*tau*d2';
+                    penalty = -obj.Hi*tau;
+                case 'd3'
+                    a = obj.alpha;
+
+                    sig = -s*a*e;
+
+                    closure = obj.Hi*sig*d3';
+                    penalty = -obj.Hi*sig;
 
                 otherwise % Unknown, boundary condition
                     error('No such boundary condition: type = %s',type);

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+scheme/TODO.txt	Tue Jul 24 20:01:29 2018 -0700
@@ -0,0 +1,1 @@
+% TODO: Rename package and abstract class to diffOp

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+scheme/bcSetup.m	Tue Jul 24 20:01:29 2018 -0700
@@ -0,0 +1,75 @@
+% function [closure, S] = bcSetup(diffOp, bc)
+% Takes a diffOp and a cell array of boundary condition definitions.
+% Each bc is a struct with the fields
+%  * type     -- Type of boundary condition
+%  * boundary -- Boundary identifier
+%  * data     -- A function_handle for a function which provides boundary data.(see below)
+% Also takes S_sign which modifies the sign of S, [-1,1]
+% Returns a closure matrix and a forcing function S.
+%
+% The boundary data function can either be a function of time or a function of time and space coordinates.
+% In the case where it only depends on time it should return the data as grid function for the boundary.
+% In the case where it also takes space coordinates the number of space coordinates should match the number of dimensions of the problem domain.
+% For example in the 2D case: f(t,x,y).
+function [closure, S] = bcSetup(diffOp, bc, S_sign)
+    default_arg('S_sign', 1);
+    assertType(bc, 'cell');
+    assert(S_sign == 1 || S_sign == -1, 'S_sign must be either 1 or -1');
+
+
+    % Setup storage arrays
+    closure = spzeros(size(diffOp));
+    gridDataPenalties = {};
+    gridDataFunctions = {};
+    symbolicDataPenalties = {};
+    symbolicDataFunctions = {};
+    symbolicDataCoords = {};
+
+    % Collect closures, penalties and data
+    for i = 1:length(bc)
+        assertType(bc{i}, 'struct');
+        [localClosure, penalty] = diffOp.boundary_condition(bc{i}.boundary, bc{i}.type);
+        closure = closure + localClosure;
+
+        if ~isfield(bc{i},'data') || isempty(bc{i}.data)
+            % Skip to next loop if there is no data
+            continue
+        end
+        assertType(bc{i}.data, 'function_handle');
+
+        % Find dimension
+        dim = size(diffOp.grid.getBoundary(bc{i}.boundary), 2);
+
+        if nargin(bc{i}.data) == 1
+            % Grid data
+            boundarySize = [size(diffOp.grid.getBoundary(bc{i}.boundary),1),1];
+            assertSize(bc{i}.data(0), boundarySize); % Eval for t = 0 and make sure the function returns a grid vector of the correct size.
+            gridDataPenalties{end+1} = penalty;
+            gridDataFunctions{end+1} = bc{i}.data;
+        elseif nargin(bc{i}.data) == 1+dim
+            % Symbolic data
+            coord = diffOp.grid.getBoundary(bc{i}.boundary);
+            symbolicDataPenalties{end+1} = penalty;
+            symbolicDataFunctions{end+1} = bc{i}.data;
+            symbolicDataCoords{end+1} = num2cell(coord ,1);
+        else
+            error('sbplib:scheme:bcSetup:DataWrongNumberOfArguments', 'bc{%d}.data has the wrong number of input arguments. Must be either only time or time and space.', i);
+        end
+    end
+
+    % Setup penalty function
+    O = spzeros(size(diffOp),1);
+    function v = S_fun(t)
+        v = O;
+        for i = 1:length(gridDataFunctions)
+            v = v + gridDataPenalties{i}*gridDataFunctions{i}(t);
+        end
+
+        for i = 1:length(symbolicDataFunctions)
+            v = v + symbolicDataPenalties{i}*symbolicDataFunctions{i}(t, symbolicDataCoords{i}{:});
+        end
+
+        v = S_sign * v;
+    end
+    S = @S_fun;
+end

--- a/+time/SBPInTimeImplicitFormulation.m	Tue Jul 24 19:59:54 2018 -0700
+++ b/+time/SBPInTimeImplicitFormulation.m	Tue Jul 24 20:01:29 2018 -0700
@@ -46,7 +46,6 @@
                 obj.f = @(t)sparse(length(v0),1);
             end
 
-
             obj.k = k;
             obj.blockSize = blockSize;
             obj.N = length(v0);
@@ -99,7 +98,7 @@
         function obj = step(obj)
             RHS = zeros(obj.blockSize*obj.N,1);
 
-            for i = 1:length(obj.blockSize)
+            for i = 1:obj.blockSize
                 RHS((1 + (i-1)*obj.N):(i*obj.N)) = obj.f(obj.t + obj.nodes(i));
             end
 

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+time/SBPInTimeScaled.m	Tue Jul 24 20:01:29 2018 -0700
@@ -0,0 +1,139 @@
+classdef SBPInTimeScaled < time.Timestepper
+    % The SBP in time method.
+    % Implemented for A*v_t = B*v + f(t), v(0) = v0
+    % The resulting system of equations is
+    %   M*u_next= K*u_prev_end + f
+    properties
+        A,B
+        f
+
+        k % total time step.
+
+        blockSize % number of points in each block
+        N % Number of components
+
+        order
+        nodes
+
+        Mtilde,Ktilde     % System matrices
+        L,U,p,q % LU factorization of M
+        e_T
+
+        scaling
+        S, Sinv % Scaling matrices
+
+        % Time state
+        t
+        vtilde
+        n
+    end
+
+    methods
+        function obj = SBPInTimeScaled(A, B, f, k, t0, v0, scaling, TYPE, order, blockSize)
+            default_arg('TYPE','gauss');
+            default_arg('f',[]);
+
+            if(strcmp(TYPE,'gauss'))
+                default_arg('order',4)
+                default_arg('blockSize',4)
+            else
+                default_arg('order', 8);
+                default_arg('blockSize',time.SBPInTimeImplicitFormulation.smallestBlockSize(order,TYPE));
+            end
+
+            obj.A = A;
+            obj.B = B;
+            obj.scaling = scaling;
+
+            if ~isempty(f)
+                obj.f = f;
+            else
+                obj.f = @(t)sparse(length(v0),1);
+            end
+
+            obj.k = k;
+            obj.blockSize = blockSize;
+            obj.N = length(v0);
+
+            obj.n = 0;
+            obj.t = t0;
+
+            %==== Build the time discretization matrix =====%
+            switch TYPE
+                case 'equidistant'
+                    ops = sbp.D2Standard(blockSize,{0,obj.k},order);
+                case 'optimal'
+                    ops = sbp.D1Nonequidistant(blockSize,{0,obj.k},order);
+                case 'minimal'
+                    ops = sbp.D1Nonequidistant(blockSize,{0,obj.k},order,'minimal');
+                case 'gauss'
+                    ops = sbp.D1Gauss(blockSize,{0,obj.k});
+            end
+
+            I = speye(size(A));
+            I_t = speye(blockSize,blockSize);
+
+            D1 = kron(ops.D1, I);
+            HI = kron(ops.HI, I);
+            e_0 = kron(ops.e_l, I);
+            e_T = kron(ops.e_r, I);
+            obj.nodes = ops.x;
+
+            % Convert to form M*w = K*v0 + f(t)
+            tau = kron(I_t, A) * e_0;
+            M = kron(I_t, A)*D1 + HI*tau*e_0' - kron(I_t, B);
+
+            K = HI*tau;
+
+            obj.S =    kron(I_t, spdiag(scaling));
+            obj.Sinv = kron(I_t, spdiag(1./scaling));
+
+            obj.Mtilde = obj.Sinv*M*obj.S;
+            obj.Ktilde = obj.Sinv*K*spdiag(scaling);
+            obj.e_T = e_T;
+
+
+            % LU factorization
+            [obj.L,obj.U,obj.p,obj.q] = lu(obj.Mtilde, 'vector');
+
+            obj.vtilde = (1./obj.scaling).*v0;
+        end
+
+        function [v,t] = getV(obj)
+            v = obj.scaling.*obj.vtilde;
+            t = obj.t;
+        end
+
+        function obj = step(obj)
+            forcing = zeros(obj.blockSize*obj.N,1);
+
+            for i = 1:obj.blockSize
+                forcing((1 + (i-1)*obj.N):(i*obj.N)) = obj.f(obj.t + obj.nodes(i));
+            end
+
+            RHS = obj.Sinv*forcing + obj.Ktilde*obj.vtilde;
+
+            y = obj.L\RHS(obj.p);
+            z = obj.U\y;
+
+            w = zeros(size(z));
+            w(obj.q) = z;
+
+            obj.vtilde = obj.e_T'*w;
+
+            obj.t = obj.t + obj.k;
+            obj.n = obj.n + 1;
+        end
+    end
+
+    methods(Static)
+        function N = smallestBlockSize(order,TYPE)
+            default_arg('TYPE','gauss')
+
+            switch TYPE
+                case 'gauss'
+                    N = 4;
+            end
+        end
+    end
+end

--- a/+time/SBPInTimeSecondOrderFormImplicit.m	Tue Jul 24 19:59:54 2018 -0700
+++ b/+time/SBPInTimeSecondOrderFormImplicit.m	Tue Jul 24 20:01:29 2018 -0700
@@ -14,15 +14,16 @@
         % Solves A*u_tt + B*u_t + C*u = f(t)
         % A, B can either both be constants or both be function handles,
         % They can also be omitted by setting them equal to the empty matrix.
-        function obj = SBPInTimeSecondOrderFormImplicit(A, B, C, f, k, t0, v0, v0t, TYPE, order, blockSize)
+        function obj = SBPInTimeSecondOrderFormImplicit(A, B, C, f, k, t0, v0, v0t, do_scaling, TYPE, order, blockSize)
             default_arg('f', []);
             default_arg('TYPE', []);
             default_arg('order', []);
             default_arg('blockSize',[]);
+            default_arg('do_scaling', false);
 
             m = length(v0);
 
-            default_arg('A', sparse(m, m));
+            default_arg('A', speye(m, m));
             default_arg('B', sparse(m, m));
             default_arg('C', sparse(m, m));
 
@@ -56,7 +57,12 @@
             obj.t = t0;
             obj.n = 0;
 
-            obj.firstOrderTimeStepper = time.SBPInTimeImplicitFormulation(obj.AA, obj.BB, obj.ff, obj.k, obj.t, w0, TYPE, order, blockSize);
+            if do_scaling
+                scaling = [ones(m,1); sqrt(diag(C))];
+                obj.firstOrderTimeStepper = time.SBPInTimeScaled(obj.AA, obj.BB, obj.ff, obj.k, obj.t, w0, scaling, TYPE, order, blockSize);
+            else
+                obj.firstOrderTimeStepper = time.SBPInTimeImplicitFormulation(obj.AA, obj.BB, obj.ff, obj.k, obj.t, w0, TYPE, order, blockSize);
+            end
         end
 
         function [v,t] = getV(obj)

--- a/+time/Timestepper.m	Tue Jul 24 19:59:54 2018 -0700
+++ b/+time/Timestepper.m	Tue Jul 24 20:01:29 2018 -0700
@@ -62,6 +62,7 @@
 
 
         function [v, t] = stepTo(obj, n, progress_bar)
+            assertScalar(n);
             default_arg('progress_bar',false);
 
             [v, t] = obj.stepN(n-obj.n, progress_bar);

--- a/.hgtags	Tue Jul 24 19:59:54 2018 -0700
+++ b/.hgtags	Tue Jul 24 20:01:29 2018 -0700
@@ -1,1 +1,2 @@
 18c023aaf3f79cbe2b9b1cf547d80babdaa1637d v0.1
+0776fa4754ff0c1918f6e1278c66f48c62d05736 grids0.1

--- a/Color.m	Tue Jul 24 19:59:54 2018 -0700
+++ b/Color.m	Tue Jul 24 20:01:29 2018 -0700
@@ -10,6 +10,10 @@
         black     = [0.000 0.000 0.000];
         white     = [1.000 1.000 1.000];
         colors = { Color.blue, Color.red, Color.yellow, Color.green, Color.purple, Color.lightblue, Color.darkred, Color.black, Color.white};
+        markers = {'+', 'o', '*', '.', 'x', 'square', 'diamond', 'v', '^', '>', '<', 'pentagram', 'hexagram'};
+        lineStyles = {'-', '--', ':', '-.'};
+
+        solidMarkers = {'o', 'square', 'diamond', 'v', 'pentagram', '^', '>', '<', 'hexagram'};
 
         notabilityYellow     = [100.0   99.0    22.0    ]/100;
         notabilityOrange     = [97.0    61.0    15.0    ]/100;
@@ -34,13 +38,11 @@
 
     methods(Static)
         function sample()
-            markers ={'+', 'o', '*', '.', 'x', 'square', 'diamond', 'v', '^', '>', '<', 'pentagram', 'hexagram'};
             % Filled and non-filled markers?
-            lineStyles = {'-', '--', ':', '-.'};
 
 
             function showMarkers(x0, y0, lx, ly, color, filled)
-                n = length(markers);
+                n = length(Color.markers);
                 s = ceil(sqrt(n));
 
                 x = linspace(x0, x0 + lx, s);
@@ -50,7 +52,7 @@
 
                 for i = 1:n
                     lh = line(X(i),Y(i));
-                    lh.Marker = markers{i};
+                    lh.Marker = Color.markers{i};
                     lh.MarkerSize = 12;
                     lh.Color = color;
 
@@ -79,13 +81,13 @@
             end
 
             function showLines(y0, ly, A, w)
-                n = length(lineStyles);
+                n = length(Color.lineStyles);
                 x = linspace(0,1,100);
                 y = linspace(y0, y0+ ly, n);
                 for i = 1:n
                     lh = line(x, y(i) + A*sin(pi*x*w));
                     lh.LineWidth = 2;
-                    lh.LineStyle = lineStyles{i};
+                    lh.LineStyle = Color.lineStyles{i};
                 end
             end
 

--- a/TextTable.m	Tue Jul 24 19:59:54 2018 -0700
+++ b/TextTable.m	Tue Jul 24 20:01:29 2018 -0700
@@ -4,28 +4,36 @@
         fmtArray
         vertDiv
         horzDiv
-
-        nCols
-        nRows
     end
 
     methods
-        function obj = TextTable(data, vertDiv, horzDiv);
+        function obj = TextTable(data, vertDiv, horzDiv)
             default_arg('vertDiv', []);
             default_arg('horzDiv', []);
 
-
             obj.data = data;
             obj.vertDiv = vertDiv;
             obj.horzDiv = horzDiv;
 
-            [obj.nRows, obj.nCols] = size(data);
             obj.fmtArray = cell(size(data));
             obj.formatAll('%s');
 
         end
 
+        function n = nRows(obj)
+            n = size(obj.data, 1);
+        end
+
+        function n = nCols(obj)
+            n = size(obj.data, 2);
+        end
+
+        function print(obj)
+            disp(obj.toString());
+        end
+
         function formatAll(obj, fmt)
+            obj.fmtArray = cell(size(obj.data));
             obj.fmtArray(:,:) = {fmt};
         end
 
@@ -58,28 +66,31 @@
 
             str = '';
 
+            N = size(strArray, 2);
+
             % First horzDiv
-            if ismember(0, obj.horzDiv)
+            if isDiv(0, obj.horzDiv, N);
                 str = [str, obj.getHorzDiv(widths)];
             end
 
             for i = 1:obj.nRows
-                str = [str, TextTable.rowToString(strArray(i,:), widths, obj.vertDiv, obj.horzDiv)];
+                str = [str, TextTable.rowToString(strArray(i,:), widths, obj.vertDiv)];
 
                 % Interior horzDiv
-                if ismember(i, obj.horzDiv)
+                if isDiv(i, obj.horzDiv, N)
                     str = [str, obj.getHorzDiv(widths)];
                 end
             end
         end
 
         function str = getHorzDiv(obj, widths)
-            str = TextTable.rowToString(cell(1,obj.nCols), widths, obj.vertDiv, obj.horzDiv);
+            str = TextTable.rowToString(cell(1,obj.nCols), widths, obj.vertDiv);
             str(find(' ' == str)) = '-';
             str(find('|' == str)) = '+';
         end
 
         function strArray = getStringArray(obj)
+            assert(all(size(obj.data) == size(obj.fmtArray)), 'Sizes of format matrix and data matrix do not match.')
             strArray = cell(size(obj.data));
 
             for i = 1:obj.nRows
@@ -91,32 +102,42 @@
     end
 
     methods (Static)
-        function str = rowToString(strs, widths, vertDiv, horzDiv)
+        function str = rowToString(strs, widths, vertDiv)
+            N = length(strs);
+
             % First vertDiv
-            if ismember(0, vertDiv)
-                str = '| ';
+            if isDiv(0, vertDiv, N)
+                prefix = '| ';
             else
-                str = ' ';
+                prefix = ' ';
             end
 
-            % Interior cols
-            for j = 1:length(strs) - 1
-                str = [str, sprintf('%*s ', widths(j), strs{j})];
+            % Pad strings
+            for i = 1:N
+                strs{i} = sprintf('%*s', widths(i), strs{i});
+            end
 
-                % Interior vertDiv
-                if ismember(j, vertDiv)
-                    str = [str, '| '];
+            % Column delimiters
+            delims = cell(1,N-1);
+            for i = 1:length(delims)
+                if isDiv(i, vertDiv, N);
+                    delims{i} = '| ';
+                else
+                    delims{i} = ' ';
                 end
             end
 
-            % Last col
-            str = [str, sprintf('%*s ', widths(end), strs{end})];
-
-            if ismember(length(strs), vertDiv)
-                str = [str, '|'];
+            if isDiv(N, vertDiv, N);
+                suffix = '|';
+            else
+                suffix = '';
             end
 
-            str = [str, sprintf('\n')];
+            str = [prefix, strjoin(strs, delims), suffix, sprintf('\n')];
         end
     end
+end
+
+function b = isDiv(i, div, N)
+    b = ismember(i, div) || ismember(i, N+div+1);
 end
\ No newline at end of file

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/assertIsMember.m	Tue Jul 24 20:01:29 2018 -0700
@@ -0,0 +1,3 @@
+function assertIsMember(v, allowed)
+    assert(ismember(v, allowed), 'Expected ''%s'' to be in the set %s', inputname(1), toString(allowed));
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/assertNumberOfArguments.m	Tue Jul 24 20:01:29 2018 -0700
@@ -0,0 +1,5 @@
+function assertNumberOfArguments(fun, N)
+    if nargin(fun) ~= N
+        error('sbplib:assertNumberOfArguments:wrongNumberOfArguments', '"%s" must have %d, found %d', inputname(1), N, nargin(fun));
+    end
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/assertScalar.m	Tue Jul 24 20:01:29 2018 -0700
@@ -0,0 +1,5 @@
+function assertScalar(obj)
+    if ~isscalar(obj)
+        error('sbplib:assertScalar:notScalar', '"%s" must be scalar, found size "%s"', inputname(1), toString(size(obj)));
+    end
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/assertSize.m	Tue Jul 24 20:01:29 2018 -0700
@@ -0,0 +1,16 @@
+% Assert that array A has the size s.
+function assertSize(A,varargin)
+    if length(varargin) == 1
+        s = varargin{1};
+        errmsg = sprintf('Expected %s to have size %s, got: %s',inputname(1), toString(s), toString(size(A)));
+        assert(all(size(A) == s), errmsg);
+    elseif length(varargin) == 2
+        dim = varargin{1};
+        s = varargin{2};
+
+        errmsg = sprintf('Expected %s to have size %d along dimension %d, got: %d',inputname(1), s, dim, size(A,dim));
+        assert(size(A,dim) == s, errmsg);
+    else
+        error('Expected 2 or 3 arguments to assertSize()');
+    end
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/assertStructFields.m	Tue Jul 24 20:01:29 2018 -0700
@@ -0,0 +1,12 @@
+% Assert that the struct s has the all the field names in the cell array fns.
+function assertStructFields(s, fns)
+    assertType(s, 'struct');
+    assertType(fns, 'cell');
+
+    ok = ismember(fns, fieldnames(s));
+    if ~all(ok)
+        str1 = sprintf("'%s' must have the fields %s\n", inputname(1), toString(fns));
+        str2 = sprintf("The following fields are missing: %s", toString(fns(~ok)));
+        error(str1 + str2);
+    end
+end

--- a/assert_size.m	Tue Jul 24 19:59:54 2018 -0700
+++ b/assert_size.m	Tue Jul 24 20:01:29 2018 -0700
@@ -1,16 +1,5 @@
 % Assert that array A has the size s.
 function assert_size(A,s)
-    errmsg = sprintf('Expected %s to have size %s, got: %s',inputname(1), format_vector(s), format_vector(size(A)));
-    assert(all(size(A) == s),errmsg);
-end
-
-function str = format_vector(a)
-    l = length(a);
-    str = sprintf('[%d',a(1));
-
-    for i = 2:l
-        str = [str sprintf(', %d',a(i))];
-    end
-
-    str = [str ']'];
+    warning('Use assertSize() instead!')
+    assertSize(A,s);
 end
\ No newline at end of file

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/convergencePlot.m	Tue Jul 24 20:01:29 2018 -0700
@@ -0,0 +1,55 @@
+function hand = convergencePlot(orders, h, e)
+    N = length(orders);
+
+    fh = figure();
+    ah = axes();
+    ah.XScale = 'log';
+    ah.YScale = 'log';
+    hold on
+    ph = {};
+    phc = {};
+    legends = {};
+    for i = 1:N
+        ph{i} = loglog(h{i}, e{i});
+        phc{i} = plotConvergenceFit(orders{i}, h{i}, e{i});
+
+        ph{i}.LineStyle = 'none';
+        ph{i}.Marker = Color.solidMarkers{i};
+        ph{i}.MarkerSize = 12;
+        ph{i}.Color = Color.colors{i};
+        ph{i}.MarkerFaceColor = Color.colors{i};
+
+        legends{i} = sprintf('$o = %d$', orders{i});
+    end
+    hold off
+
+    lh = legend([ph{:}], legends);
+    lh.Interpreter = 'latex';
+    lh.Location = 'SouthEast';
+
+    for i = 1:N
+        uistack(phc{i}, 'bottom');
+    end
+
+    xlabel('$h$', 'interpreter', 'latex')
+    ylabel('Error', 'interpreter', 'latex')
+
+    % xlim([0.7e-2, 1e-1])
+    % ylim([3e-5, 4])
+
+    grid on
+
+    ah = gca();
+    ah.TickLabelInterpreter = 'latex';
+    setFontSize(fh);
+
+    % if savePngs
+    %     savepng(fh, 'fig/conv/conv',600)
+    % end
+
+    hand = struct();
+    hand.fig = fh;
+    hand.data = ph;
+    hand.fits = phc;
+    hand.legend = lh;
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/four.m	Tue Jul 24 20:01:29 2018 -0700
@@ -0,0 +1,19 @@
+% four returns the fourier transform u_hat of the function u and the frequencies w
+function [w, u_hat] = four(x, u)
+    u_hat = fft(u);
+
+    N = length(x);
+    L = x(end) - x(1);
+
+    k = shift_k(0:N-1);
+
+    u_hat = fftshift(u_hat);
+
+    dw = 2*pi/L;
+    w = dw*k;
+end
+
+function k_shifted = shift_k(k)
+    N = length(k);
+    k_shifted = [-floor(N/2):-1, 0, 1:ceil(N/2)-1];
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/fourInv.m	Tue Jul 24 20:01:29 2018 -0700
@@ -0,0 +1,4 @@
+function u = ifour(u_hat)
+    u_hat = ifftshift(u_hat);
+    u = ifft(u_hat);
+end

--- a/gaussian.m	Tue Jul 24 19:59:54 2018 -0700
+++ b/gaussian.m	Tue Jul 24 20:01:29 2018 -0700
@@ -1,3 +1,3 @@
 function z = gaussian(x,x0,d)
-    z = exp(-norm(x-x0).^2/d^2);
+    z = exp(-sum((x-x0).^2,2)/d^2);
 end
\ No newline at end of file

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/matlabFunctionSizePreserving.m	Tue Jul 24 20:01:29 2018 -0700
@@ -0,0 +1,12 @@
+% Takes a symfun and makes a better anonymous function
+function fun = matlabFunctionSizePreserving(f)
+    mf = matlabFunction(f);
+    args = argnames(f);
+
+    funStr = func2str(mf);
+    for i = 1:length(args)
+        funStr = [funStr sprintf(' + 0*%s', toString(args(i)))];
+    end
+
+    fun = str2func(funStr);
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mononomial.m	Tue Jul 24 20:01:29 2018 -0700
@@ -0,0 +1,17 @@
+% calculate a N-D mononomial with powers k in points x:
+%  z = x(:,1).^k(1) * x(:,2).^k(2) * ...
+function z = mononomial(x, k)
+    assert(size(x,2) == length(k), 'k must have the same length as the width of x');
+
+    if any(k < 0)
+        z = x(:,1)*0;
+        return
+    end
+
+    denom = prod(factorial(k));
+
+    for i = 1:length(k)
+        x(:,i) = x(:,i).^k(i);
+    end
+    z = prod(x,2)/denom;
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/nextColor.m	Tue Jul 24 20:01:29 2018 -0700
@@ -0,0 +1,5 @@
+function c = nextColor(ah)
+    default_arg('ah', gca);
+
+    c = ah.ColorOrder(ah.ColorOrderIndex, :);
+end
\ No newline at end of file

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/pointIndex.m	Tue Jul 24 20:01:29 2018 -0700
@@ -0,0 +1,4 @@
+% Get the index of the points p within the tall array of points ps
+function [I, ok] = pointIndex(p, ps)
+    [ok, I] = ismember(p,  ps, 'rows');
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/rickerWavelet.m	Tue Jul 24 20:01:29 2018 -0700
@@ -0,0 +1,3 @@
+function y = rickerWavelet(x, x0, A)
+    y = (1-2*pi^2*A^2*(x-x0).^2).*exp(-pi^2*A^2*(x-x0).^2);
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/sbplibLocation.m	Tue Jul 24 20:01:29 2018 -0700
@@ -0,0 +1,4 @@
+function location = sbplibLocation()
+    scriptname  = mfilename('fullpath');
+    [location, ~, ~] = fileparts(scriptname);
+end

--- a/sbplibVersion.m	Tue Jul 24 19:59:54 2018 -0700
+++ b/sbplibVersion.m	Tue Jul 24 20:01:29 2018 -0700
@@ -1,11 +1,10 @@
 % Prints the version and location of the sbplib currently in use.
 function sbplibVersion()
-    scriptname  = mfilename('fullpath');
-    [folder,~,~] = fileparts(scriptname);
+    location = sbplibLocation();
 
     name = 'sbplib (feature/grids)';
     ver = '0.0.x';
 
     fprintf('%s %s\n', name, ver);
-    fprintf('Running in:\n%s\n',folder);
+    fprintf('Running in:\n%s\n', location);
 end
\ No newline at end of file

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/stencilEquation.m	Tue Jul 24 20:01:29 2018 -0700
@@ -0,0 +1,13 @@
+% Find the equation for the stencil for d^k/dx^k
+function [A,b] = stencilEquation(k, offsets, order)
+    q = sym('q', [1, length(offsets)]);
+
+    p = 0:(order-1+k);
+
+    v     = vandermonde(offsets, p);
+    vdiff = vandermonde(      0, p-k);
+
+    eq = q*v == vdiff;
+
+    [A,b] = equationsToMatrix(eq, q);
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/subsSymfun.m	Tue Jul 24 20:01:29 2018 -0700
@@ -0,0 +1,14 @@
+% Subs for a symfun
+% f remains a symbolic function. If any of it's arguments is eliminated
+% it is removed from the argument list while preserving the order of the
+% other arguments
+function f = subsSymfun(f, old, new)
+    args = argnames(f);
+
+    newExpr = subs(f, old, new);
+    vars = symvar(subs(args, old, new));
+
+    newArgs = args(ismember(args,vars));
+
+    f = symfun(newExpr, newArgs);
+end
\ No newline at end of file

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/vandermonde.m	Tue Jul 24 20:01:29 2018 -0700
@@ -0,0 +1,15 @@
+% Create vandermonde matrix for points x and polynomials of order p
+% x is a list of N points of size [N,dim],
+% p is a list of polynomial orders of size [M, dim].
+% the given mononomials are evaluated and the NxM matrix V is returned.
+function V = vandermonde(x, p)
+    assert(size(x,2) == size(p,2), 'x and p must have the same number of columns')
+    n = size(x,1);
+    m = size(p,1);
+
+    for i = 1:m
+        V(:,i) = mononomial(x, p(i,:));
+    end
+
+    assertSize(V,[n,m]);
+end