diff --git a/ext/MTKBifurcationKitExt.jl b/ext/MTKBifurcationKitExt.jl
index b709b2eec0..6bc536958f 100644
--- a/ext/MTKBifurcationKitExt.jl
+++ b/ext/MTKBifurcationKitExt.jl
@@ -144,9 +144,10 @@ function BifurcationKit.BifurcationProblem(osys::ODESystem, args...; kwargs...)
     if !ModelingToolkit.iscomplete(osys)
         error("A completed `ODESystem` is required. Call `complete` or `structural_simplify` on the system before creating a `BifurcationProblem`")
     end
-    nsys = NonlinearSystem([0 ~ eq.rhs for eq in equations(osys)],
+    nsys = NonlinearSystem([0 ~ eq.rhs for eq in full_equations(osys)],
         unknowns(osys),
         parameters(osys);
+        observed = observed(osys),
         name = nameof(osys))
     return BifurcationKit.BifurcationProblem(complete(nsys), args...; kwargs...)
 end
diff --git a/test/extensions/bifurcationkit.jl b/test/extensions/bifurcationkit.jl
index 861d8f45f9..698dd085c8 100644
--- a/test/extensions/bifurcationkit.jl
+++ b/test/extensions/bifurcationkit.jl
@@ -133,3 +133,33 @@ let
     @test length(fold_points) == 2
     @test fold_points ≈ [-1.1851851706940317, -5.6734983580551894e-6] # test that they occur at the correct parameter values).
 end
+
+let
+    @mtkmodel FOL begin
+        @parameters begin
+            τ # parameters
+        end
+        @variables begin
+            x(t) # dependent variables
+            RHS(t)
+        end
+        @equations begin
+            RHS ~ τ + x^2 - 0.1
+            D(x) ~ RHS
+        end
+    end
+
+    @mtkbuild fol = FOL()
+
+    par = [fol.τ => 0.0]
+    u0 = [fol.x => -1.0]
+    #prob = ODEProblem(fol, u0, (0.0, 1.), par)
+
+    bif_par = fol.τ
+    bp = BifurcationProblem(fol, u0, par, bif_par)
+    opts_br = ContinuationPar(p_min = -1.0,
+        p_max = 1.0)
+    bf = bifurcationdiagram(bp, PALC(), 2, opts_br)
+
+    @test bf.γ.specialpoint[1].param≈0.1 atol=1e-4 rtol=1e-4
+end