--- a/+draw/prompt_point.m	Sun Jun 17 14:44:05 2018 -0700
+++ b/+draw/prompt_point.m	Sun Jun 17 15:29:46 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	Sun Jun 17 14:44:05 2018 -0700
+++ b/+grid/Cartesian.m	Sun Jun 17 15:29:46 2018 -0700
@@ -16,7 +16,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/+multiblock/Laplace.m	Sun Jun 17 15:29:46 2018 -0700
@@ -0,0 +1,52 @@
+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 * obj.jacobian();
+        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 [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

--- a/+scheme/bcSetup.m	Sun Jun 17 14:44:05 2018 -0700
+++ b/+scheme/bcSetup.m	Sun Jun 17 15:29:46 2018 -0700
@@ -3,43 +3,70 @@
 % Each bc is a struct with the fields
 %  * type     -- Type of boundary condition
 %  * boundary -- Boundary identifier
-%  * data     -- A function_handle with time and space coordinates as a parameters, for example f(t,x,y) for a 2D problem
+%  * 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
+% 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));
-    penalties = {};
-    dataFunctions = {};
-    dataParams = {};
+    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');
 
-        coord = diffOp.grid.getBoundary(bc{i}.boundary);
-        assertNumberOfArguments(bc{i}.data, 1+size(coord,2));
+        % Find dimension
+        dim = size(diffOp.grid.getBoundary(bc{i}.boundary), 2);
 
-        penalties{end+1} = penalty;
-        dataFunctions{end+1} = bc{i}.data;
-        dataParams{end+1} = num2cell(coord ,1);
+        if nargin(bc{i}.data) == 1
+            % Grid data
+            boundarySize = [size(diffOp.grid.getBoundary(bc{i}.boundary),1),1];
+            assert_size(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(dataFunctions)
-            v = v + penalties{i}*dataFunctions{i}(t, dataParams{i}{:});
+        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;

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+time/SBPInTimeScaled.m	Sun Jun 17 15:29:46 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	Sun Jun 17 14:44:05 2018 -0700
+++ b/+time/SBPInTimeSecondOrderFormImplicit.m	Sun Jun 17 15:29:46 2018 -0700
@@ -14,11 +14,12 @@
         % 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', true);
 
             m = length(v0);
 
@@ -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)