diff --git a/src/Blocks/Blocks.jl b/src/Blocks/Blocks.jl
index 08bc163e7..5610bbb8f 100644
--- a/src/Blocks/Blocks.jl
+++ b/src/Blocks/Blocks.jl
@@ -16,7 +16,7 @@ export Abs, Sign, Sqrt, Sin, Cos, Tan, Asin, Acos, Atan, Atan2, Sinh, Cosh, Tanh
 export Log, Log10
 include("math.jl")
 
-export Constant, Sine, Cosine, ContinuousClock, Ramp, Step, ExpSine
+export Constant, Sine, Cosine, ContinuousClock, Ramp, Step, ExpSine, Square, Triangular
 include("sources.jl")
 
 export Limiter, DeadZone, SlewRateLimiter
diff --git a/src/Blocks/sources.jl b/src/Blocks/sources.jl
index 3daa331e5..315cb0ff5 100644
--- a/src/Blocks/sources.jl
+++ b/src/Blocks/sources.jl
@@ -1,3 +1,47 @@
+# Define and register smooth functions
+# These are "smooth" aka differentiable and avoid Gibbs effect
+# These follow: `offset` + `smooth_wave` * `smooth_step` with zero output for `t < start_time`
+function smooth_cos(x, δ, f, amplitude, ϕ, offset, start_time)
+    offset + amplitude * cos(2*π*f*(x - start_time) + ϕ) * smooth_step(x, δ, one(x), zero(x), start_time)
+end
+
+function smooth_damped_sin(x, δ, f, amplitude, damping, ϕ, offset, start_time)
+    offset + exp((start_time - x)*damping)*amplitude*sin(2*π*f*(x - start_time) + ϕ) * smooth_step(x, δ, one(x), zero(x), start_time)
+end
+
+function smooth_ramp(x, δ, height, duration, offset, start_time)
+    offset + height/(duration) * (smooth_xH(x, δ, start_time) - smooth_xH(x, δ, start_time+duration))
+end
+
+function smooth_sin(x, δ, f, amplitude, ϕ, offset, start_time)
+    offset + amplitude * sin(2*pi*f*(x - start_time) + ϕ) * smooth_step(x, δ, one(x), zero(x), start_time)
+end
+
+function smooth_square(x, δ, f, amplitude, offset, start_time)
+    offset + amplitude*2atan(sin(2π*(x - start_time)*f)/δ)/π * smooth_step(x, δ, one(x), zero(x), start_time)
+end
+
+function smooth_step(x, δ, height, offset, start_time)
+    offset + height*(atan((x - start_time)/δ)/π + 0.5)
+end
+
+function smooth_triangular(x, δ, f, amplitude, offset, start_time)
+    offset + amplitude * (1-2acos((1 - δ)sin(2π*(x - start_time)*f))/π) * smooth_step(x, δ, one(x), zero(x), start_time)
+end
+
+function smooth_xH(x, δ, tₒ)
+    0.5*(x-tₒ) * (1+((x-tₒ)/sqrt((x-tₒ)^2+δ^2)))
+end
+
+function square(x, f, amplitude, offset, start_time)
+    offset + (x > start_time) * (amplitude * (4*floor(f*(x - start_time)) - 2*floor(2*(x - start_time)*f) + 1))
+end
+
+function triangular(x, f, amplitude, offset, start_time)
+    p = 1/f # period
+    offset + (x > start_time) * (4 * amplitude * f * abs(abs((x - p/4 - start_time) % p) - p/2) - amplitude)
+end
+
 """
 Generate constant signal.
 
@@ -22,25 +66,36 @@ Generate sine signal.
 # Parameters:
 - `frequency`: [Hz] Frequency of sine wave
 - `amplitude`: Amplitude of sine wave
-- `phase`: [rad] Phase of sine wave 
+- `phase`: [rad] Phase of sine wave
 - `offset`: Offset of output signal
 - `start_time`: [s] Output `y = offset` for `t < start_time`
+- `smooth`:  If `true`, returns a smooth wave. Defaults to `false`
+             It uses a smoothing factor of `δ=1e-5`
 
 # Connectors:
 - `output`
 """
-function Sine(;name, 
-    frequency, 
+function Sine(;name,
+    frequency,
     amplitude=1,
     phase=0,
     offset=0,
-    start_time=0)
+    start_time=0,
+    smooth=false)
 
     @named output = RealOutput()
     pars = @parameters offset=offset start_time=start_time amplitude=amplitude frequency=frequency phase=phase
+    equation = if smooth == false
+        offset + ifelse(t < start_time, 0, amplitude* sin(2*pi*frequency*(t - start_time) + phase))
+    else
+        δ = 1e-5
+        smooth_sin(t, δ, frequency, amplitude, phase, offset, start_time)
+    end
+
     eqs = [
-        output.u ~ offset + ifelse(t < start_time, 0, amplitude* sin(2*pi*frequency*(t - start_time) + phase))
+        output.u ~ equation
     ]
+
     compose(ODESystem(eqs, t, [], pars; name=name), [output])
 end
 
@@ -50,32 +105,43 @@ Generate cosine signal.
 # Parameters:
 - `frequency`: [Hz] Frequency of sine wave
 - `amplitude`: Amplitude of sine wave
-- `phase`: [rad] Phase of sine wave 
+- `phase`: [rad] Phase of sine wave
 - `offset`: Offset of output signal
 - `start_time`: [s] Output `y = offset` for `t < start_time`
+- `smooth`:  If `true`, returns a smooth wave. Defaults to `false`
+             It uses a smoothing factor of `δ=1e-5`
 
 # Connectors:
 - `output`
 """
-function Cosine(;name, 
-    frequency, 
+
+function Cosine(;name,
+    frequency,
     amplitude=1,
     phase=0,
     offset=0,
-    start_time=0)
+    start_time=0,
+    smooth=false)
 
     @named output = RealOutput()
     pars = @parameters offset=offset start_time=start_time amplitude=amplitude frequency=frequency phase=phase
+    equation = if smooth == false
+        offset + ifelse(t < start_time, zero(t), amplitude* cos(2*pi*frequency*(t - start_time) + phase))
+    else
+        δ = 1e-5
+        smooth_cos(t, δ, frequency, amplitude, phase, offset, start_time)
+    end
     eqs = [
-        output.u ~ offset + ifelse(t < start_time, 0, amplitude* cos(2*pi*frequency*(t - start_time) + phase))
+        output.u ~ equation
     ]
+
     compose(ODESystem(eqs, t, [], pars; name=name), [output])
 end
 
 """
 Generate current time signal.
 
-# Parameters: 
+# Parameters:
 - `offset`: Offset of output signal
 - `start_time`: [s] Output `y = offset` for `t < start_time`
 
@@ -86,8 +152,9 @@ function ContinuousClock(;name, offset=0, start_time=0)
     @named output = RealOutput()
     pars = @parameters offset=offset start_time=start_time
     eqs = [
-        output.u ~ offset + ifelse(t < start_time, 0, t - start_time)
+        output.u ~ offset + ifelse(t < start_time, zero(t), t - start_time)
     ]
+
     compose(ODESystem(eqs, t, [], pars; name=name), [output])
 end
 
@@ -99,22 +166,73 @@ Generate ramp signal.
 - `duration`: [s] Duration of ramp (= 0.0 gives a Step)
 - `offset`: Offset of output signal
 - `start_time`: [s] Output `y = offset` for `t < start_time`
+- `smooth`:  If `true`, returns a smooth wave. Defaults to `false`
+             It uses a smoothing factor of `δ=1e-5`
 
 # Connectors:
 - `output`
 """
-function Ramp(;name, 
-    offset=0,
+function Ramp(;name,
     height=1,
-    duration=1, 
-    start_time=0)
+    duration=1,
+    offset=0,
+    start_time=0,
+    smooth=false)
 
     @named output = RealOutput()
     pars = @parameters offset=offset start_time=start_time height=height duration=duration
-    eqs = [
-        output.u ~ offset + ifelse(t < start_time, 0, 
+    equation = if smooth == false
+        offset + ifelse(t < start_time, 0,
             ifelse(t < (start_time + duration), (t - start_time) * height / duration, height))
+    else
+        δ = 1e-5
+        smooth_ramp(t, δ, height, duration, offset, start_time)
+    end
+
+    eqs = [
+        output.u ~ equation
+    ]
+
+    compose(ODESystem(eqs, t, [], pars; name=name), [output])
+end
+
+"""
+Generate smooth square signal.
+
+# Parameters:
+- `frequency`: [Hz] Frequency of square wave
+- `amplitude`: Amplitude of square wave
+- `offset`: Offset of output signal
+- `start_time`: [s] Output `y = offset` for `t < start_time`
+- `smooth`:  If `true`, returns a smooth wave. Defaults to `false`
+             It uses a smoothing factor of `δ=1e-5`
+
+# Connectors:
+- `output`
+"""
+function Square(; name, frequency=1.0, amplitude=1.0,
+    offset=0.0, start_time=0.0, smooth=false)
+    δ = 1e-5
+
+    @named output = RealOutput()
+    pars = @parameters begin
+        frequency=frequency
+        amplitude=amplitude
+        offset=offset
+        start_time=start_time
+    end
+
+    equation = if smooth == false
+        square(t, frequency, amplitude, offset, start_time)
+    else
+        δ = 1e-5
+        smooth_square(t, δ, frequency, amplitude, offset, start_time)
+    end
+
+    eqs = [
+        output.u ~ equation
     ]
+
     compose(ODESystem(eqs, t, [], pars; name=name), [output])
 end
 
@@ -125,16 +243,26 @@ Generate step signal.
 - `height`: Height of step
 - `offset`: Offset of output signal
 - `start_time`: [s] Output `y = offset` for `t < start_time`
+- `smooth`:  If `true`, returns a smooth wave. Defaults to `false`
+             It uses a smoothing factor of `δ=1e-5`
 
 # Connectors:
 - `output`
 """
-function Step(;name, offset=0, height=1, start_time=0)
+function Step(;name, height=1, offset=0, start_time=0, smooth=true)
     @named output = RealOutput()
     pars = @parameters offset=offset start_time=start_time height=height
+    equation = if smooth == false
+        offset + ifelse(t < start_time, zero(t), height)
+    else
+        δ = 1e-5
+        smooth_step(t, δ, height, offset, start_time)
+    end
+
     eqs = [
-        output.u ~ offset + ifelse(t < start_time, 0, height)
+        output.u ~ equation
     ]
+
     compose(ODESystem(eqs, t, [], pars; name=name), [output])
 end
 
@@ -145,26 +273,77 @@ Generate exponentially damped sine signal.
 - `frequency`: [Hz] Frequency of sine wave
 - `amplitude`: Amplitude of sine wave
 - `damping`: [1/s] Damping coefficient of sine wave
-- `phase`: [rad] Phase of sine wave 
+- `phase`: [rad] Phase of sine wave
 - `offset`: Offset of output signal
 - `start_time`: [s] Output `y = offset` for `t < start_time`
+- `smooth`:  If `true`, returns a smooth wave. Defaults to `false`
+             It uses a smoothing factor of `δ=1e-5`
 
 # Connectors:
 - `output`
 """
-function ExpSine(;name, 
-    frequency, 
+function ExpSine(; name,
+    frequency,
     amplitude=1,
     damping=0.1,
     phase=0,
     offset=0,
-    start_time=0)
+    start_time=0,
+    smooth=false)
 
     @named output = RealOutput()
     pars = @parameters offset=offset start_time=start_time amplitude=amplitude frequency=frequency phase=phase damping=damping
+
+    equation = if smooth == false
+        offset + ifelse(t < start_time, 0, amplitude * exp(-damping * (t - start_time)) * sin(2*pi*frequency*(t - start_time) + phase))
+    else
+        δ = 1e-5
+        smooth_damped_sin(t, δ, frequency, amplitude, damping, phase, offset, start_time)
+    end
+
     eqs = [
-        output.u ~ offset + ifelse(t < start_time, 0, amplitude * exp(-damping * (t - start_time)) * sin(2*pi*frequency*(t - start_time) + phase))
+        output.u ~ equation
     ]
+
+    compose(ODESystem(eqs, t, [], pars; name=name), [output])
+end
+
+"""
+Generate smooth triangular signal for frequencies less than or equal to 25 Hz
+
+# Parameters:
+- `frequency`: [Hz] Frequency of square wave
+- `amplitude`: Amplitude of square wave
+- `offset`: Offset of output signal.
+- `start_time`: [s] Output `y = offset` for `t < start_time`
+- `smooth`:  If `true`, returns a smooth wave. Defaults to `false`
+             It uses a smoothing factor of `δ=1e-5`
+
+# Connectors:
+- `output`
+"""
+function Triangular(; name, amplitude=1.0, frequency=1.0,
+    offset=0.0, start_time=0.0, smooth=false)
+
+    @named output = RealOutput()
+    pars = @parameters begin
+        amplitude=amplitude
+        frequency=frequency
+        offset=offset
+        start_time=start_time
+    end
+
+    equation = if smooth == false
+        triangular(t, frequency, amplitude, offset, start_time)
+    else
+        δ = 1e-5
+        smooth_triangular(t, δ, frequency, amplitude, offset, start_time)
+    end
+
+    eqs = [
+        output.u ~ equation
+    ]
+
     compose(ODESystem(eqs, t, [], pars; name=name), [output])
 end
 
diff --git a/test/Blocks/sources.jl b/test/Blocks/sources.jl
index 2ac7f22cf..16a21ee04 100644
--- a/test/Blocks/sources.jl
+++ b/test/Blocks/sources.jl
@@ -1,5 +1,6 @@
 using ModelingToolkit, ModelingToolkitStandardLibrary, OrdinaryDiffEq
 using ModelingToolkitStandardLibrary.Blocks
+using ModelingToolkitStandardLibrary.Blocks: smooth_sin, smooth_cos, smooth_damped_sin, smooth_square, smooth_step, smooth_ramp, smooth_triangular, triangular, square
 
 @parameters t
 
@@ -27,10 +28,11 @@ end
     amplitude=2
     phase=0
     offset=1
-    start_time=0
+    start_time=2
+    δ = 1e-5
+    @named int = Integrator()
 
     @named src = Sine(frequency=frequency, amplitude=amplitude, phase=phase, offset=offset, start_time=start_time)
-    @named int = Integrator()
     @named iosys = ODESystem([
         connect(src.output, int.input),
         ],
@@ -43,36 +45,81 @@ end
 
     sol = solve(prob, Rodas4())
     @test sol[src.output.u] ≈ sine.(sol.t, frequency, amplitude, phase, offset, start_time) atol=1e-3
+
+    @named smooth_src = Sine(frequency=frequency,
+        amplitude=amplitude,
+        phase=phase,
+        offset=offset,
+        start_time=start_time,
+        smooth=true)
+    @named smooth_iosys = ODESystem([
+        connect(smooth_src.output, int.input),
+        ],
+        t,
+        systems=[int, smooth_src],
+    )
+
+    smooth_sys = structural_simplify(smooth_iosys)
+    smooth_prob = ODEProblem(smooth_sys, Pair[int.x=>0.0], (0.0, 10.0))
+    smooth_sol = solve(smooth_prob, Rodas4())
+
+    @test smooth_sol[smooth_src.output.u] ≈ smooth_sin.(smooth_sol.t, δ, frequency, amplitude, phase, offset, start_time) atol=1e-3
 end
 
 @testset "Cosine" begin
-    cosine(t, frequency, amplitude, phase, offset, start_time) = offset + ifelse(t < start_time, 0, amplitude* cos(2*pi*frequency*(t - start_time) + phase))
-    
+    cosine(t, frequency, amplitude, phase, offset, start_time) =
+        offset + ifelse(t < start_time, 0, amplitude* cos(2*pi*frequency*(t - start_time) + phase))
+
     frequency=1
     amplitude=2
     phase=0
     offset=1
-    start_time=0
-
-    @named src = Cosine(frequency=frequency, amplitude=amplitude, phase=phase, offset=offset, start_time=start_time)
+    start_time=2
+    δ = 1e-5
     @named int = Integrator()
+
+    @named src = Cosine(frequency=frequency,
+        amplitude=amplitude,
+        phase=phase,
+        offset=offset,
+        start_time=start_time,
+        smooth=false)
     @named iosys = ODESystem([
         connect(src.output, int.input),
         ],
         t,
         systems=[int, src],
     )
-    sys = structural_simplify(iosys)
 
+    sys = structural_simplify(iosys)
     prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 10.0))
-
     sol = solve(prob, Rodas4())
+
     @test sol[src.output.u] ≈ cosine.(sol.t, frequency, amplitude, phase, offset, start_time) atol=1e-3
+
+    @named smooth_src = Cosine(frequency=frequency,
+        amplitude=amplitude,
+        phase=phase,
+        offset=offset,
+        start_time=start_time,
+        smooth=true)
+    @named smooth_iosys = ODESystem([
+        connect(smooth_src.output, int.input),
+        ],
+        t,
+        systems=[int, smooth_src],
+    )
+
+    smooth_sys = structural_simplify(smooth_iosys)
+    smooth_prob = ODEProblem(smooth_sys, Pair[int.x=>0.0], (0.0, 10.0))
+    smooth_sol = solve(smooth_prob, Rodas4())
+
+    @test smooth_sol[smooth_src.output.u] ≈ smooth_cos.(smooth_sol.t, δ, frequency, amplitude, phase, offset, start_time) atol=1e-3
 end
 
 @testset "ContinuousClock" begin
     cont_clock(t, offset, start_time) = offset + ifelse(t < start_time, 0, t - start_time)
-    
+
     offset, start_time = 1, 0
 
     @named src = ContinuousClock(offset=offset, start_time=start_time)
@@ -93,51 +140,158 @@ end
 
 @testset "Ramp" begin
     ramp(t, offset, height, duration, start_time) = offset + ifelse(t < start_time, 0, ifelse(t < (start_time + duration), (t - start_time) * height / duration, height))
-    
-    offset, height, duration, start_time = 1, 2, 2, 0
 
-    @named src = Ramp(offset=offset, height=height, duration=duration, start_time=start_time)
+    offset, height, duration, start_time, δ = 1, 2, 2, 0, 1e-5
     @named int = Integrator()
+
+    @named src = Ramp(offset=offset, height=height, duration=duration, start_time=start_time)
     @named iosys = ODESystem([
         connect(src.output, int.input),
         ],
         t,
         systems=[int, src],
     )
-    sys = structural_simplify(iosys)
 
+    sys = structural_simplify(iosys)
     prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 10.0))
-
     sol = solve(prob, Rodas4())
+
     @test sol[src.output.u] ≈ ramp.(sol.t, offset, height, duration, start_time) atol=1e-3
+
+    start_time = 2
+    @named smooth_src = Ramp(offset=offset, height=height, duration=duration, start_time=start_time, smooth=true)
+    @named smooth_iosys = ODESystem([
+        connect(smooth_src.output, int.input),
+        ],
+        t,
+        systems=[int, smooth_src],
+    )
+
+    smooth_sys = structural_simplify(smooth_iosys)
+    smooth_prob = ODEProblem(smooth_sys, Pair[int.x=>0.0], (0.0, 10.0))
+    smooth_sol = solve(smooth_prob, Rodas4())
+
+    @test smooth_sol[smooth_src.output.u] ≈ smooth_ramp.(smooth_sol.t, δ, height, duration, offset, start_time) atol=1e-3
 end
 
 @testset "Step" begin
     step(t, offset, height, start_time) = offset + ifelse(t < start_time, 0, height)
 
-    offset, height, start_time = 1, 2, 5
+    offset, height, start_time, δ = 1, 2, 5, 1e-5
+    @named int = Integrator()
 
     @named src = Step(offset=offset, height=height, start_time=start_time)
-    @named int = Integrator()
     @named iosys = ODESystem([
         connect(src.output, int.input),
         ],
         t,
         systems=[int, src],
     )
-    sys = structural_simplify(iosys)
 
+    sys = structural_simplify(iosys)
     prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 10.0))
+    sol = solve(prob, Rodas4())
+
+    @test sol[src.output.u] ≈ step.(sol.t, offset, height, start_time) atol=1e-2
+
+    @named smooth_src = Step(offset=offset, height=height, start_time=start_time, smooth=true)
+    @named smooth_iosys = ODESystem([
+        connect(smooth_src.output, int.input),
+        ],
+        t,
+        systems=[int, smooth_src],
+    )
+
+    smooth_sys = structural_simplify(smooth_iosys)
+    smooth_prob = ODEProblem(smooth_sys, Pair[int.x=>0.0], (0.0, 10.0))
+    smooth_sol = solve(smooth_prob, Rodas4())
+
+    @test smooth_sol[smooth_src.output.u] ≈ smooth_step.(smooth_sol.t, δ, height, offset, start_time)
+end
+
+@testset "Square" begin
+    frequency=1
+    amplitude=2
+    offset=1
+    start_time=2.5
+    δ = 1e-5
+    @named int = Integrator()
 
+    @named src = Square(frequency=frequency, amplitude=amplitude,
+        offset=offset, start_time=start_time)
+    @named iosys = ODESystem([
+        connect(src.output, int.input),
+        ],
+        t,
+        systems=[int, src],
+    )
+
+    sys = structural_simplify(iosys)
+    prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 10.0))
     sol = solve(prob, Rodas4())
-    @test sol[src.output.u] ≈ step.(sol.t, offset, height, start_time) atol=1e-3
+
+    @test sol[src.output.u] ≈ square.(sol.t, frequency, amplitude, offset, start_time) atol=1e-3
+
+    @named smooth_src = Square(frequency=frequency, amplitude=amplitude,
+    offset=offset, start_time=start_time, smooth=true)
+    @named smooth_iosys = ODESystem([
+        connect(smooth_src.output, int.input),
+        ],
+        t,
+        systems=[int, smooth_src]
+    )
+
+    smooth_sys = structural_simplify(smooth_iosys)
+    smooth_prob = ODEProblem(smooth_sys, Pair[int.x=>0.0], (0.0, 10.0))
+    smooth_sol = solve(smooth_prob, Rodas4())
+
+    @test smooth_sol[smooth_src.output.u] ≈ smooth_square.(smooth_sol.t, δ, frequency, amplitude, offset, start_time) atol=1e-3
+end
+
+@testset "Triangular" begin
+    frequency=5
+    amplitude=1
+    offset = 2
+    start_time = 1
+    δ = 1e-5
+    @named int = Integrator()
+
+    @named src = Triangular(frequency=frequency, amplitude=amplitude,
+        offset=offset, start_time=start_time)
+    @named iosys = ODESystem([
+        connect(src.output, int.input),
+        ],
+        t,
+        systems=[int, src],
+    )
+
+    sys = structural_simplify(iosys)
+    prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 4.0))
+    sol = solve(prob, Rodas4(), saveat=0.01)
+
+    @test sol[src.output.u] ≈ triangular.(sol.t, frequency, amplitude, offset, start_time) atol=1e-3
+
+    @named smooth_src = Triangular(frequency=frequency, amplitude=amplitude,
+        offset=offset, start_time=start_time, smooth=true)
+    @named smooth_iosys = ODESystem([
+        connect(smooth_src.output, int.input),
+        ],
+        t,
+        systems=[int, smooth_src],
+    )
+
+    smooth_sys = structural_simplify(smooth_iosys)
+    smooth_prob = ODEProblem(smooth_sys, Pair[int.x=>0.0], (0.0, 4.0))
+    smooth_sol = solve(smooth_prob, Rodas4(), saveat=0.01)
+
+    @test smooth_sol[smooth_src.output.u] ≈ smooth_triangular.(smooth_sol.t, δ, frequency, amplitude, offset, start_time) atol=1e-3
 end
 
 @testset "ExpSine" begin
     exp_sine(t, amplitude, frequency, damping, phase, start_time) = offset + ifelse(t < start_time, 0, amplitude * exp(-damping * (t - start_time)) * sin(2*pi*frequency*(t - start_time) + phase))
 
-    frequency, amplitude, damping, phase, offset, start_time = 3, 2, 0.10, 0, 0, 0
-
+    frequency, amplitude, damping = 3, 2, 0.10
+    phase, offset, start_time, δ =  0, 0, 0, 1e-5
     @named src = ExpSine(frequency=frequency, amplitude=amplitude, damping=damping, phase=phase, offset=offset, start_time=start_time)
     @named int = Integrator()
     @named iosys = ODESystem([
@@ -147,9 +301,22 @@ end
         systems=[int, src],
     )
     sys = structural_simplify(iosys)
-
     prob = ODEProblem(sys, Pair[int.x=>0.0], (0.0, 10.0))
-
     sol = solve(prob, Rodas4())
+
     @test sol[src.output.u] ≈ exp_sine.(sol.t, amplitude, frequency, damping, phase, start_time) atol=1e-3
-end
\ No newline at end of file
+
+    offset, start_time = 1, 2
+    @named smooth_src = ExpSine(frequency=frequency, amplitude=amplitude, damping=damping, phase=phase, offset=offset, start_time=start_time, smooth=true)
+    @named smooth_iosys = ODESystem([
+        connect(smooth_src.output, int.input),
+        ],
+        t,
+        systems=[int, smooth_src],
+    )
+    smooth_sys = structural_simplify(smooth_iosys)
+    smooth_prob = ODEProblem(smooth_sys, Pair[int.x=>0.0], (0.0, 10.0))
+    smooth_sol = solve(smooth_prob, Rodas4())
+
+    @test smooth_sol[smooth_src.output.u] ≈ smooth_damped_sin.(smooth_sol.t, δ, frequency, amplitude, damping, phase, offset, start_time) atol=1e-3
+end
diff --git a/test/Thermal/thermal.jl b/test/Thermal/thermal.jl
index 59b739b01..794493fd9 100644
--- a/test/Thermal/thermal.jl
+++ b/test/Thermal/thermal.jl
@@ -32,7 +32,7 @@ D = Differential(t)
 
     # Check if Relative temperature sensor reads the temperature of heat capacitor
     # when connected to a thermal conductor and a fixed temperature source
-    @test_broken sol[reltem_sensor.T] + sol[tem_src.port.T] == sol[mass1.T] + sol[th_conductor.dT]
+    @test sol[reltem_sensor.T] + sol[tem_src.port.T] == sol[mass1.T] + sol[th_conductor.dT]
 
     @info "Building a two-body system..."
     eqs = [
@@ -58,8 +58,8 @@ D = Differential(t)
     m1, m2 = sol.u[end]
     @test m1 ≈ m2 atol=1e-1
     mass_T = reduce(hcat, sol.u)
-    @test_broken sol[T_sensor1.T] == mass_T[1, :]
-    @test_broken sol[T_sensor2.T] == mass_T[2, :]
+    @test sol[T_sensor1.T] == mass_T[1, :]
+    @test sol[T_sensor2.T] == mass_T[2, :]
 end
 
 # Test HeatFlowSensor, FixedHeatFlow, ThermalResistor, ThermalConductor
@@ -92,11 +92,11 @@ end
     prob = DAEProblem(sys, D.(states(sys)) .=> 0.0, u0, (0, 3.0))
     sol  = solve(prob, DFBDF())
 
-    @test_broken sol[th_conductor.dT] .* G == sol[th_conductor.Q_flow]
-    @test_broken sol[th_conductor.Q_flow] ≈ sol[hf_sensor1.Q_flow] + sol[flow_src.port.Q_flow]
+    @test sol[th_conductor.dT] .* G == sol[th_conductor.Q_flow]
+    @test sol[th_conductor.Q_flow] ≈ sol[hf_sensor1.Q_flow] + sol[flow_src.port.Q_flow]
 
-    @test_broken sol[mass1.T] == sol[th_resistor.port_a.T]
-    @test_broken sol[th_resistor.dT] ./ R ≈ sol[th_resistor.Q_flow]
+    @test sol[mass1.T] == sol[th_resistor.port_a.T]
+    @test sol[th_resistor.dT] ./ R ≈ sol[th_resistor.Q_flow]
 
 end
 
@@ -130,7 +130,7 @@ end
 
     # Heat-flow-rate is equal in magnitude
     # and opposite in direction
-    @test_broken sol[gas.Q_flow] + sol[coolant.Q_flow] == zeros(length(sol))
+    @test sol[gas.Q_flow] + sol[coolant.Q_flow] == zeros(length(sol))
 end
 
 # Test ConvectiveConductor, BodyRadiation
@@ -163,9 +163,9 @@ end
     prob = DAEProblem(sys, D.(states(sys)) .=> 0.0, u0, (0, 3.0))
     sol = solve(prob, DFBDF())
 
-    @test_broken sol[dissipator.dT] == sol[radiator.port_a.T] - sol[radiator.port_b.T]
+    @test sol[dissipator.dT] == sol[radiator.port_a.T] - sol[radiator.port_b.T]
     rad_Q_flow = G*σ*(Tᵧ^4 - Tᵪ^4)
-    @test_broken sol[radiator.Q_flow] == fill(rad_Q_flow, length(sol[radiator.Q_flow]))
+    @test sol[radiator.Q_flow] == fill(rad_Q_flow, length(sol[radiator.Q_flow]))
 end
 
 @testset "Thermal Collector" begin
@@ -197,7 +197,7 @@ end
     prob = DAEProblem(sys, D.(states(sys)) .=> 0.0, u0, (0, 3.0))
     sol  = solve(prob, DFBDF())
 
-    @test_broken sol[collector.port_b.Q_flow] + sol[collector.port_a1.Q_flow] + sol[collector.port_a2.Q_flow] ==
+    @test sol[collector.port_b.Q_flow] + sol[collector.port_a1.Q_flow] + sol[collector.port_a2.Q_flow] ==
         zeros(length(sol[collector.port_b.Q_flow]))
-    @test_broken sol[collector.port_b.T] == sol[collector.port_a1.T] == sol[collector.port_a2.T]
+    @test sol[collector.port_b.T] == sol[collector.port_a1.T] == sol[collector.port_a2.T]
 end