--- 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)