diff -r 7ab7d42a5b24 -r 3a286d2d1939 +multiblock/+domain/Annulus.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+multiblock/+domain/Annulus.m	Mon Feb 14 11:14:46 2022 +0100
@@ -0,0 +1,74 @@
+classdef Annulus < multiblock.DefCurvilinear
+    properties
+        r_inner % Radii of inner disk
+        c_inner % Center of inner disk
+        r_outer % Radii of outer disk
+        c_outer % Radii of outer disk
+    end
+
+    methods
+        function obj = Annulus(r_inner, c_inner, r_outer, c_outer)
+            default_arg('r_inner', 0.3);
+            default_arg('c_inner', [0; 0]);
+            default_arg('r_outer', 1)
+            default_arg('c_outer', [0; 0]);
+            % Assert that the problem is well-defined
+            d = norm(c_outer-c_inner,2);
+            assert(r_inner > 0, 'Inner radius must be greater than zero');
+            assert(r_outer > d+r_inner, 'Inner disk not contained in outer disk');
+            
+            cir_out_A = parametrization.Curve.circle(c_outer,r_outer,[-pi/2 pi/2]);
+            cir_in_A = parametrization.Curve.circle(c_inner,r_inner,[pi/2 -pi/2]);
+            
+            cir_out_B = parametrization.Curve.circle(c_outer,r_outer,[pi/2 3*pi/2]);
+            cir_in_B = parametrization.Curve.circle(c_inner,r_inner,[3*pi/2 pi/2]);
+
+            c0_out = cir_out_A(0);
+            c1_out = cir_out_A(1);
+            
+            c0_in_A = cir_in_A(1);
+            c1_in_A = cir_in_A(0);
+            
+            c0_out_B = cir_out_B(0);
+            c1_out_B = cir_out_B(1);
+            
+            c0_in_B = cir_in_B(1);
+            c1_in_B = cir_in_B(0);
+
+
+            sp2_A = parametrization.Curve.line(c0_in_A,c0_out);
+            sp3_A = parametrization.Curve.line(c1_in_A,c1_out);
+            
+            sp2_B = parametrization.Curve.line(c0_in_B,c0_out_B);
+            sp3_B = parametrization.Curve.line(c1_in_B,c1_out_B);
+
+
+            A = parametrization.Ti(sp2_A, cir_out_A, sp3_A.reverse, cir_in_A); 
+            B = parametrization.Ti(sp2_B , cir_out_B,sp3_B.reverse, cir_in_B );
+            
+            blocks = {A,B};
+            blocksNames = {'A','B'};
+
+            conn = cell(2,2);
+
+            conn{1,2} = {'n','s'};
+            conn{2,1} = {'n','s'};
+
+            boundaryGroups = struct();
+            boundaryGroups.out = multiblock.BoundaryGroup({{1,'e'},{2,'e'}});
+            boundaryGroups.in = multiblock.BoundaryGroup({{1,'w'},{2,'w'}});
+            boundaryGroups.all = multiblock.BoundaryGroup({{1,'e'},{2,'w'},{1,'w'},{2,'e'}});
+
+            obj = obj@multiblock.DefCurvilinear(blocks, conn, boundaryGroups, blocksNames);
+
+            obj.r_inner = r_inner;
+            obj.r_outer = r_outer;
+            obj.c_inner = c_inner;
+            obj.c_outer = c_outer;
+        end
+
+        function ms = getGridSizes(obj, m)
+            ms = {[m m], [m m]};
+        end
+    end
+end

diff -r 7ab7d42a5b24 -r 3a286d2d1939 +time/Rungekutta4SecondOrder.m
--- a/+time/Rungekutta4SecondOrder.m	Mon Feb 14 10:49:49 2022 +0100
+++ b/+time/Rungekutta4SecondOrder.m	Mon Feb 14 11:14:46 2022 +0100
@@ -36,9 +36,32 @@
             obj.S = S;
             obj.m = length(v0);
             obj.n = 0;
-
+            
+            % Construct RHS system on first order form. 
+            % 3 different setups are handeled:
+            % If none of D, E, S are time-dependent (e.g function_handles)
+            % then a complete matrix of the RHS is constructed.
+            % If the source term S is time-depdentent, then the first order system matrix M
+            % is formed via D and E, while S is kept as a function handle
+            % If D or E are time-dependent, then all terns are treated as function handles.
+            if ~isa(D, 'function_handle') && ~isa(E, 'function_handle') && isa(S, 'function_handle')
+                default_arg('D', sparse(obj.m, obj.m));
+                default_arg('E', sparse(obj.m, obj.m));
+                default_arg('S', @(t)sparse(obj.m, 1));
+                
+                I = speye(obj.m);
+                O = sparse(obj.m,obj.m);
 
-            if isa(D, 'function_handle') || isa(E, 'function_handle') || isa(S, 'function_handle')
+                obj.M = [
+                    O, I;
+                    D, E;
+                ];
+
+                obj.k = k;
+                obj.t = t0;
+                obj.w = [v0; v0t];
+                obj.F = @(w,t)obj.rhsTimeDepS(w,t,obj.M,S,obj.m);
+            elseif isa(D, 'function_handle') || isa(E, 'function_handle')
                 default_arg('D', @(t)sparse(obj.m, obj.m));
                 default_arg('E', @(t)sparse(obj.m, obj.m));
                 default_arg('S', @(t)sparse(obj.m, 1)    );
@@ -63,7 +86,6 @@
                     D(t)*w(1:obj.m) + E(t)*w(obj.m+1:end) + S(t);
                 ];
             else
-
                 default_arg('D', sparse(obj.m, obj.m));
                 default_arg('E', sparse(obj.m, obj.m));
                 default_arg('S', sparse(obj.m, 1)    );
@@ -110,6 +132,11 @@
         function k = getTimeStep(lambda)
             k = rk4.get_rk4_time_step(lambda);
         end
+        
+        function F = rhsTimeDepS(w,t,M,S,m)
+            F = M*w;
+            F(m+1:end) = F(m+1:end) + S(t);
+        end
     end
 
 end
\ No newline at end of file