diff --git a/docs/pages.jl b/docs/pages.jl
index 866134eaf..4746259a8 100644
--- a/docs/pages.jl
+++ b/docs/pages.jl
@@ -4,6 +4,7 @@ pages = [
         "RC Circuit" => "tutorials/rc_circuit.md",
         "Custom Components" => "tutorials/custom_component.md",
         "Thermal Conduction Model" => "tutorials/thermal_model.md",
+        "DC Motor with Speed Controller" => "tutorials/dc_motor_pi.md",
     ],
     "API" => [
         "Basic Blocks" => "API/blocks.md",
diff --git a/docs/src/index.md b/docs/src/index.md
index 5c39fa999..a1f0417c0 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -17,6 +17,7 @@ import Pkg; Pkg.add("ModelingToolkitStandardLibrary")
 - [RC Circuit](http://mtkstdlib.sciml.ai/dev/tutorials/rc_circuit/)
 - [Custom Component](http://mtkstdlib.sciml.ai/dev/tutorials/custom_component/)
 - [Thermal Model](http://mtkstdlib.sciml.ai/dev/tutorials/thermal_model/)
+- [DC Motor with PI-controller](http://mtkstdlib.sciml.ai/dev/tutorials/dc_motor_pi/)
 
 ## Libraries
 
diff --git a/docs/src/tutorials/dc_motor_pi.md b/docs/src/tutorials/dc_motor_pi.md
new file mode 100644
index 000000000..41ccd077a
--- /dev/null
+++ b/docs/src/tutorials/dc_motor_pi.md
@@ -0,0 +1,93 @@
+# DC Motor with PI-controller
+In this example a PI-controller is setup for speed control of a DC-motor. First the needed packages 
+are imported and the parameters of the model defined.
+
+```@example dc_motor_pi
+using ModelingToolkit
+using ModelingToolkitStandardLibrary.Electrical
+using ModelingToolkitStandardLibrary.Mechanical.Rotational
+using ModelingToolkitStandardLibrary.Blocks
+using OrdinaryDiffEq
+using Plots
+
+@parameters t
+
+R = 0.5 # [Ohm]
+L = 4.5e-3 # [H]
+k = 0.5 # [N.m/A]
+J = 0.02 # [kg.m²]
+f = 0.01 # [N.m.s/rad]
+tau_L_step = -3 # [N.m]
+nothing # hide
+```
+
+The actual model can now be composed.
+
+```@example dc_motor_pi
+@named ground = Ground()
+@named source = Voltage()
+@named ref = Blocks.Step(height = 1, start_time = 0)
+@named pi_controller = Blocks.LimPI(k = 1.1, T = 0.035, u_max = 0.6, Ta = 0.035)
+@named feedback = Blocks.Feedback()
+@named R1 = Resistor(R = R)
+@named L1 = Inductor(L = L)
+@named emf = EMF(k = k)
+@named fixed = Fixed()
+@named load = Torque(use_support = false)
+@named load_step = Blocks.Step(height = tau_L_step, start_time = 3)
+@named inertia = Inertia(J = J)
+@named friction = Damper(d = f)
+@named speed_sensor = SpeedSensor()
+
+connections = [connect(fixed.flange, emf.support, friction.flange_b)
+               connect(emf.flange, friction.flange_a, inertia.flange_a)
+               connect(inertia.flange_b, load.flange)
+               connect(inertia.flange_b, speed_sensor.flange)
+               connect(load_step.output, load.tau)
+               connect(ref.output, feedback.input1)
+               connect(speed_sensor.w, feedback.input2)
+               connect(feedback.output, pi_controller.err_input)
+               connect(pi_controller.ctr_output, source.V)
+               connect(source.p, R1.p)
+               connect(R1.n, L1.p)
+               connect(L1.n, emf.p)
+               connect(emf.n, source.n, ground.g)]
+
+@named model = ODESystem(connections, t,
+                         systems = [
+                             ground,
+                             ref,
+                             pi_controller,
+                             feedback,
+                             source,
+                             R1,
+                             L1,
+                             emf,
+                             fixed,
+                             load,
+                             load_step,
+                             inertia,
+                             friction,
+                             speed_sensor,
+                         ])
+nothing # hide
+```
+
+Now the model can be simulated. Typical rotational mechanical systems are described via `DAE` 
+(differential algebraic equations), however in this case ModelingToolkit can simplify the model enough
+so that it can be represented as a system of `ODEs` (ordinary differential equations).
+
+```@example dc_motor_pi
+sys = structural_simplify(model)
+prob = ODEProblem(sys, [], (0, 6.0))
+sol = solve(prob, Rodas4())
+
+p1 = Plots.plot(sol.t, sol[inertia.w], ylabel = "Angular Vel. in rad/s",
+                label = "Measurement", title = "DC Motor with Speed Controller")
+Plots.plot!(sol.t, sol[ref.output.u], label = "Reference")
+p2 = Plots.plot(sol.t, sol[load.tau.u], ylabel = "Disturbance in Nm", label = "")
+Plots.plot(p1, p2, layout = (2, 1))
+nothing # hide
+```
+
+![DC Motor with speed controller](dc_motor_pi.png)
diff --git a/src/Electrical/Analog/ideal_components.jl b/src/Electrical/Analog/ideal_components.jl
index e4ed71698..959e3d8b7 100644
--- a/src/Electrical/Analog/ideal_components.jl
+++ b/src/Electrical/Analog/ideal_components.jl
@@ -26,7 +26,7 @@ See [OnePort](@ref)
 - `n` Negative pin
 
 # Parameters:
-- `R`: [`Ω`] Resistance
+- `R`: [`Ohm`] Resistance
 """
 function Resistor(; name, R)
     @named oneport = OnePort()
@@ -70,8 +70,7 @@ end
 Creates an ideal capacitor.
 
 # States:
-- `v(t)`: [`V`]
-  The voltage across the capacitor, given by `D(v) ~ p.i / C`
+- `v(t)`: [`V`] The voltage across the capacitor, given by `D(v) ~ p.i / C`
 
 # Connectors:
 - `p` Positive pin
@@ -160,3 +159,73 @@ function Short(; name)
     eqs = [v ~ 0]
     extend(ODESystem(eqs, t, [], []; name = name), oneport)
 end
+
+"""
+    HeatingResistor(;name, R_ref=1.0, T_ref=300.15, alpha=0)
+
+Temperature dependent electrical resistor
+
+# States
+- See [OnePort](@ref)
+- `R(t)`: [`Ohm`] Temperature dependent resistance `R ~ R_ref*(1 + alpha*(heat_port.T(t) - T_ref))`
+
+# Connectors
+- `p` Positive pin
+- `n` Negative pin
+
+# Parameters: 
+- `R_ref`: [`Ω`] Reference resistance 
+- `T_ref`: [K] Reference temperature
+"""
+function HeatingResistor(; name, R_ref = 1.0, T_ref = 300.15, alpha = 0)
+    @named oneport = OnePort()
+    @unpack v, i = oneport
+    @named heat_port = HeatPort()
+    pars = @parameters begin
+        R_ref = R_ref
+        T_ref = T_ref
+        alpha = alpha
+    end
+    @variables R(t) = R_ref
+    eqs = [R ~ R_ref * (1 + alpha * (heat_port.T - T_ref))
+           heat_port.Q_flow ~ -v * i # -LossPower
+           v ~ i * R]
+    extend(ODESystem(eqs, t, [R], pars; name = name, systems = [heat_port]), oneport)
+end
+
+"""
+    EMF(;name, k)
+
+Electromotoric force (electric/mechanic transformer)
+
+# States
+- `v(t)`: [`V`] The voltage across component `p.v - n.v`
+- `i(t)`: [`A`] The current passing through positive pin
+- `phi`: [`rad`] Rotation angle (=flange.phi - support.phi)
+- `w`: [`rad/s`] Angular velocity (= der(phi))
+
+# Connectors
+- `p` [Pin](@ref) Positive pin
+- `n` [Pin](@ref) Negative pin
+- `flange` [Flange](@ref) Shaft of EMF shaft
+- `support` [Support](@ref) Support/housing of emf shaft
+
+# Parameters: 
+- `k`: [`N⋅m/A`] Transformation coefficient 
+"""
+function EMF(; name, k)
+    @named p = Pin()
+    @named n = Pin()
+    @named flange = Flange()
+    @named support = Support()
+    @parameters k = k
+    @variables v(t)=0.0 i(t)=0.0 phi(t)=0.0 w(t)=0.0
+    eqs = [v ~ p.v - n.v
+           0 ~ p.i + n.i
+           i ~ p.i
+           phi ~ flange.phi - support.phi
+           D(phi) ~ w
+           k * w ~ v
+           flange.tau ~ -k * i]
+    ODESystem(eqs, t, [v, i, phi, w], [k]; name = name, systems = [p, n, flange, support])
+end
diff --git a/src/Electrical/Electrical.jl b/src/Electrical/Electrical.jl
index 3a394252a..544efce78 100644
--- a/src/Electrical/Electrical.jl
+++ b/src/Electrical/Electrical.jl
@@ -6,16 +6,26 @@ module Electrical
 
 using ModelingToolkit, Symbolics, IfElse
 using OffsetArrays
+using ..Thermal: HeatPort
+using ..Mechanical.Rotational: Flange, Support
+using ..Blocks: RealInput, RealOutput
 
 @parameters t
 D = Differential(t)
 
-using ..Blocks: RealInput, RealOutput
-
+export Pin, OnePort
 include("utils.jl")
+
+export Capacitor, Ground, Inductor, Resistor, Conductor, Short, IdealOpAmp, EMF,
+       HeatingResistor
 include("Analog/ideal_components.jl")
+
+export CurrentSensor, PotentialSensor, VoltageSensor, PowerSensor, MultiSensor
 include("Analog/sensors.jl")
+
+export Voltage, Current
 include("Analog/sources.jl")
+
 # include("Digital/components.jl")
 # include("Digital/gates.jl")
 # include("Digital/tables.jl")
@@ -26,22 +36,4 @@ include("Analog/sources.jl")
 # - machines
 # - multi-phase
 
-export #Interface
-       Pin,
-# Analog Components
-       Capacitor, Ground, Inductor, Resistor, Conductor,
-       Short, IdealOpAmp,
-# Analog Sensors
-       CurrentSensor, PotentialSensor, VoltageSensor,
-       PowerSensor, MultiSensor,
-# Analog Sources
-       Voltage, Current
-
-# # Digital Gates
-# And, Or, Not, Xor, Nand, Nor, Xnor,
-# # Digital components
-# HalfAdder, FullAdder, MUX, DEMUX, Encoder, Decoder,
-# # Digital Sources
-# DigitalPin, Pulse, PulseDiff
-
 end
diff --git a/src/Mechanical/Rotational/components.jl b/src/Mechanical/Rotational/components.jl
index c5cd8e050..b1c821634 100644
--- a/src/Mechanical/Rotational/components.jl
+++ b/src/Mechanical/Rotational/components.jl
@@ -39,6 +39,7 @@ end
 function Inertia(; name, J, phi_start = 0.0, w_start = 0.0, a_start = 0.0)
     @named flange_a = Flange()
     @named flange_b = Flange()
+    J > 0 || throw(ArgumentError("Expected `J` to be positive"))
     @parameters J = J
     sts = @variables begin
         phi(t) = phi_start
@@ -73,6 +74,7 @@ Linear 1D rotational spring
 function Spring(; name, c, phi_rel0 = 0.0)
     @named partial_comp = PartialCompliant()
     @unpack phi_rel, tau = partial_comp
+    c > 0 || throw(ArgumentError("Expected `c` to be positive"))
     pars = @parameters begin
         c = c
         phi_rel0 = phi_rel0
@@ -102,6 +104,7 @@ Linear 1D rotational damper
 function Damper(; name, d)
     @named partial_comp = PartialCompliantWithRelativeStates()
     @unpack w_rel, tau = partial_comp
+    d > 0 || throw(ArgumentError("Expected `d` to be positive"))
     pars = @parameters d = d
     eqs = [tau ~ d * w_rel]
     extend(ODESystem(eqs, t, [], pars; name = name), partial_comp)
@@ -130,6 +133,7 @@ This element characterizes any type of gear box which is fixed in the ground and
 function IdealGear(; name, ratio, use_support = false)
     @named partial_element = PartialElementaryTwoFlangesAndSupport2(use_support = use_support)
     @unpack phi_support, flange_a, flange_b = partial_element
+    ratio > 0 || throw(ArgumentError("Expected `ratio` to be positive"))
     @parameters ratio = ratio
     sts = @variables phi_a(t)=0.0 phi_b(t)=0.0
     eqs = [phi_a ~ flange_a.phi - phi_support
diff --git a/src/Mechanical/Rotational/sources.jl b/src/Mechanical/Rotational/sources.jl
index bd46b914f..06a375a9b 100644
--- a/src/Mechanical/Rotational/sources.jl
+++ b/src/Mechanical/Rotational/sources.jl
@@ -1,5 +1,5 @@
 """
-    Torque(;name) 
+    Torque(; name, use_support=false) 
 
 Input signal acting as external torque on a flange
 
diff --git a/src/Mechanical/Rotational/utils.jl b/src/Mechanical/Rotational/utils.jl
index 78386d87d..75ba614a6 100644
--- a/src/Mechanical/Rotational/utils.jl
+++ b/src/Mechanical/Rotational/utils.jl
@@ -118,7 +118,7 @@ Partial model for a component with one rotational 1-dim. shaft flange and a supp
 function PartialElementaryOneFlangeAndSupport2(; name, use_support = false)
     @named flange = Flange()
     sys = [flange]
-    @variables phi_support(t)
+    @variables phi_support(t) = 0.0
     if use_support
         @named support = Support()
         eqs = [support.phi ~ phi_support
diff --git a/src/ModelingToolkitStandardLibrary.jl b/src/ModelingToolkitStandardLibrary.jl
index e8a24088d..ebb910eff 100644
--- a/src/ModelingToolkitStandardLibrary.jl
+++ b/src/ModelingToolkitStandardLibrary.jl
@@ -2,8 +2,8 @@ module ModelingToolkitStandardLibrary
 
 include("Blocks/Blocks.jl")
 include("Mechanical/Mechanical.jl")
+include("Thermal/Thermal.jl")
 include("Electrical/Electrical.jl")
 include("Magnetic/Magnetic.jl")
-include("Thermal/Thermal.jl")
 
 end
diff --git a/test/Blocks/sources.jl b/test/Blocks/sources.jl
index fd53fccf7..258637ee8 100644
--- a/test/Blocks/sources.jl
+++ b/test/Blocks/sources.jl
@@ -102,7 +102,6 @@ end
     sys = structural_simplify(iosys)
     prob = ODEProblem(sys, Pair[int.x => 0.0], (0.0, 10.0))
     sol = solve(prob, Rodas4())
-
     @test sol.retcode == :Success
     @test sol[src.output.u]≈cosine.(sol.t, frequency, amplitude, phase, offset, start_time) atol=1e-3
 
@@ -169,7 +168,6 @@ end
     sys = structural_simplify(iosys)
     prob = ODEProblem(sys, Pair[int.x => 0.0], (0.0, 10.0))
     sol = solve(prob, Rodas4())
-
     @test sol.retcode == :Success
     @test sol[src.output.u]≈ramp.(sol.t, offset, height, duration, start_time) atol=1e-3
 
@@ -365,7 +363,6 @@ end
     sys = structural_simplify(iosys)
     prob = ODEProblem(sys, Pair[int.x => 0.0], (0.0, 10.0))
     sol = solve(prob, Rodas4())
-
     @test sol.retcode == :Success
     @test sol[src.output.u]≈exp_sine.(sol.t, amplitude, frequency, damping, phase,
                                       start_time) atol=1e-3
diff --git a/test/Electrical/analog.jl b/test/Electrical/analog.jl
index 74f5eab35..5e3b807a5 100644
--- a/test/Electrical/analog.jl
+++ b/test/Electrical/analog.jl
@@ -157,10 +157,10 @@ end
                                  systems = [resistor, capacitor, source, ground, voltage])
         sys = structural_simplify(model)
         prob = ODAEProblem(sys, [capacitor.v => 10.0], (0.0, 10.0))
-        sol = solve(prob, Rodas5())
-        @test sol.retcode == :Success
         sol = solve(prob, Tsit5())
         @test sol.retcode == :Success
+        sol = solve(prob, Rodas4())
+        @test sol.retcode == :Success
 
         # Plots.plot(sol; vars=[voltage.v, capacitor.v])
     end
diff --git a/test/multi_domain.jl b/test/multi_domain.jl
new file mode 100644
index 000000000..78688f408
--- /dev/null
+++ b/test/multi_domain.jl
@@ -0,0 +1,215 @@
+using ModelingToolkitStandardLibrary.Electrical
+using ModelingToolkitStandardLibrary.Mechanical.Rotational
+using ModelingToolkitStandardLibrary.Blocks
+using ModelingToolkitStandardLibrary.Blocks
+using ModelingToolkitStandardLibrary.Thermal
+import ModelingToolkitStandardLibrary
+using ModelingToolkit, OrdinaryDiffEq, Test
+# using Plots
+
+@parameters t
+D = Differential(t)
+
+@testset "DC motor" begin
+    R = 0.5
+    L = 4.5e-3
+    k = 0.5
+    J = 0.02
+    f = 0.01
+    V_step = 10
+    tau_L_step = -3
+    @named ground = Ground()
+    @named source = Voltage()
+    @named voltage_step = Blocks.Step(height = V_step, start_time = 0)
+    @named R1 = Resistor(R = R)
+    @named L1 = Inductor(L = L)
+    @named emf = EMF(k = k)
+    @named fixed = Fixed()
+    @named load = Torque(use_support = false)
+    @named load_step = Blocks.Step(height = tau_L_step, start_time = 3)
+    @named inertia = Inertia(J = J)
+    @named friction = Damper(d = f)
+
+    connections = [connect(fixed.flange, emf.support, friction.flange_b)
+                   connect(emf.flange, friction.flange_a, inertia.flange_a)
+                   connect(inertia.flange_b, load.flange)
+                   connect(load_step.output, load.tau)
+                   connect(voltage_step.output, source.V)
+                   connect(source.p, R1.p)
+                   connect(R1.n, L1.p)
+                   connect(L1.n, emf.p)
+                   connect(emf.n, source.n, ground.g)]
+
+    @named model = ODESystem(connections, t,
+                             systems = [
+                                 ground,
+                                 voltage_step,
+                                 source,
+                                 R1,
+                                 L1,
+                                 emf,
+                                 fixed,
+                                 load,
+                                 load_step,
+                                 inertia,
+                                 friction,
+                             ])
+    sys = structural_simplify(model)
+
+    @test_broken prob = ODAEProblem(sys, Pair[], (0, 6.0))
+    @test_skip begin
+        sol = solve(prob, Rodas4())
+        @test sol.retcode == :Success
+        # EMF equations
+        @test -0.5 .* sol[emf.i] == sol[emf.flange.tau]
+        @test sol[emf.v] == 0.5 .* sol[emf.w]
+        # test steady-state values
+        dc_gain = [f/(k^2 + f * R) k/(k^2 + f * R); k/(k^2 + f * R) -R/(k^2 + f * R)]
+        idx_t = findfirst(sol.t .> 2.5)
+        @test sol[inertia.w][idx_t]≈(dc_gain * [V_step; 0])[2] rtol=1e-3
+        @test sol[emf.i][idx_t]≈(dc_gain * [V_step; 0])[1] rtol=1e-3
+        idx_t = findfirst(sol.t .> 5.5)
+        @test sol[inertia.w][idx_t]≈(dc_gain * [V_step; -tau_L_step])[2] rtol=1e-3
+        @test sol[emf.i][idx_t]≈(dc_gain * [V_step; -tau_L_step])[1] rtol=1e-3
+    end
+
+    prob = DAEProblem(sys, D.(states(sys)) .=> 0.0, Pair[], (0, 6.0))
+    sol = solve(prob, DFBDF())
+    @test sol.retcode == :Success
+    # EMF equations
+    @test -0.5 .* sol[emf.i] == sol[emf.flange.tau]
+    @test sol[emf.v] == 0.5 .* sol[emf.w]
+    # test steady-state values
+    dc_gain = [f/(k^2 + f * R) k/(k^2 + f * R); k/(k^2 + f * R) -R/(k^2 + f * R)]
+    idx_t = findfirst(sol.t .> 2.5)
+    @test sol[inertia.w][idx_t]≈(dc_gain * [V_step; 0])[2] rtol=1e-3
+    @test sol[emf.i][idx_t]≈(dc_gain * [V_step; 0])[1] rtol=1e-3
+    idx_t = findfirst(sol.t .> 5.5)
+    @test sol[inertia.w][idx_t]≈(dc_gain * [V_step; -tau_L_step])[2] rtol=1e-3
+    @test sol[emf.i][idx_t]≈(dc_gain * [V_step; -tau_L_step])[1] rtol=1e-3
+
+    # p1 = Plots.plot(sol, vars=[inertia.w], ylabel="Angular Vel. in rad/s", label="")
+    # p2 = Plots.plot(sol, vars=[emf.i], ylabel="Current in A", label="")
+    # Plots.plot(p1, p2, layout=(2,1))
+    # Plots.savefig("dc_motor.png")
+end
+
+@testset "DC motor with speed sensor" begin
+    R = 0.5
+    L = 4.5e-3
+    k = 0.5
+    J = 0.02
+    f = 0.01
+    V_step = 10
+    tau_L_step = -3
+    @named ground = Ground()
+    @named source = Voltage()
+    @named voltage_step = Blocks.Step(height = V_step, start_time = 0)
+    @named R1 = Resistor(R = R)
+    @named L1 = Inductor(L = L)
+    @named emf = EMF(k = k)
+    @named fixed = Fixed()
+    @named load = Torque(use_support = false)
+    @named load_step = Blocks.Step(height = tau_L_step, start_time = 3)
+    @named inertia = Inertia(J = J)
+    @named friction = Damper(d = f)
+    @named speed_sensor = SpeedSensor()
+
+    connections = [connect(fixed.flange, emf.support, friction.flange_b)
+                   connect(emf.flange, friction.flange_a, inertia.flange_a)
+                   connect(inertia.flange_b, load.flange)
+                   connect(inertia.flange_b, speed_sensor.flange)
+                   connect(load_step.output, load.tau)
+                   connect(voltage_step.output, source.V)
+                   connect(source.p, R1.p)
+                   connect(R1.n, L1.p)
+                   connect(L1.n, emf.p)
+                   connect(emf.n, source.n, ground.g)]
+
+    @named model = ODESystem(connections, t,
+                             systems = [
+                                 ground,
+                                 voltage_step,
+                                 source,
+                                 R1,
+                                 L1,
+                                 emf,
+                                 fixed,
+                                 load,
+                                 load_step,
+                                 inertia,
+                                 friction,
+                                 speed_sensor,
+                             ])
+    sys = structural_simplify(model)
+
+    @test_broken prob = ODAEProblem(sys, Pair[], (0, 6.0)) # KeyError: key 17 not found
+    @test_skip begin
+        sol = solve(prob, Rodas4())
+
+        @test sol.retcode == :Success
+        # EMF equations
+        @test -0.5 .* sol[emf.i] == sol[emf.flange.tau]
+        @test sol[emf.v] == 0.5 .* sol[emf.w]
+
+        # test steady-state values
+        dc_gain = [f/(k^2 + f * R) k/(k^2 + f * R); k/(k^2 + f * R) -R/(k^2 + f * R)]
+        idx_t = findfirst(sol.t .> 2.5)
+        @test sol[inertia.w][idx_t]≈(dc_gain * [V_step; 0])[2] rtol=1e-3
+        @test sol[emf.i][idx_t]≈(dc_gain * [V_step; 0])[1] rtol=1e-3
+        idx_t = findfirst(sol.t .> 5.5)
+        @test sol[inertia.w][idx_t]≈(dc_gain * [V_step; -tau_L_step])[2] rtol=1e-3
+        @test sol[emf.i][idx_t]≈(dc_gain * [V_step; -tau_L_step])[1] rtol=1e-3
+
+        # 
+        @test all(sol[inertia.w] .== sol[speed_sensor.w.u])
+    end
+
+    prob = DAEProblem(sys, D.(states(sys)) .=> 0.0, Pair[], (0, 6.0))
+    sol = solve(prob, DFBDF())
+
+    @test sol.retcode == :Success
+    # EMF equations
+    @test -0.5 .* sol[emf.i] == sol[emf.flange.tau]
+    @test sol[emf.v] == 0.5 .* sol[emf.w]
+    # test steady-state values
+    dc_gain = [f/(k^2 + f * R) k/(k^2 + f * R); k/(k^2 + f * R) -R/(k^2 + f * R)]
+    idx_t = findfirst(sol.t .> 2.5)
+    @test sol[inertia.w][idx_t]≈(dc_gain * [V_step; 0])[2] rtol=1e-3
+    @test sol[emf.i][idx_t]≈(dc_gain * [V_step; 0])[1] rtol=1e-3
+    idx_t = findfirst(sol.t .> 5.5)
+    @test sol[inertia.w][idx_t]≈(dc_gain * [V_step; -tau_L_step])[2] rtol=1e-3
+    @test sol[emf.i][idx_t]≈(dc_gain * [V_step; -tau_L_step])[1] rtol=1e-3
+    # 
+    @test all(sol[inertia.w] .== sol[speed_sensor.w.u])
+end
+
+@testset "El. Heating Circuit" begin
+    @named ground = Ground()
+    @named source = Voltage()
+    @named voltage_sine = Blocks.Sine(amplitude = 220, frequency = 1)
+    @named heating_resistor = HeatingResistor(R_ref = 100, alpha = 1e-3, T_ref = 293.15)
+    @named thermal_conductor = ThermalConductor(G = 50)
+    @named env = FixedTemperature(T = 273.15 + 20)
+    connections = [connect(source.n, ground.g, heating_resistor.n)
+                   connect(source.p, heating_resistor.p)
+                   connect(voltage_sine.output, source.V)
+                   connect(heating_resistor.heat_port, thermal_conductor.port_a)
+                   connect(thermal_conductor.port_b, env.port)]
+
+    @named model = ODESystem(connections, t,
+                             systems = [
+                                 ground,
+                                 voltage_sine,
+                                 source,
+                                 heating_resistor,
+                                 thermal_conductor,
+                                 env,
+                             ])
+    sys = structural_simplify(model)
+
+    prob = ODEProblem(sys, Pair[], (0, 6.0))
+    sol = solve(prob, Rodas4())
+    @test sol.retcode == :Success
+    @test sol[source.v * source.i] == -sol[env.port.Q_flow]
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 9cf6bcbea..982465f84 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -22,3 +22,5 @@ using SafeTestsets
 
 # Mechanical
 @safetestset "Mechanical" begin include("Mechanical/rotational.jl") end
+
+@safetestset "Multi-Domain" begin include("multi_domain.jl") end