diff --git a/README.md b/README.md
index abecb955b..8ad7d97b0 100644
--- a/README.md
+++ b/README.md
@@ -11,7 +11,7 @@
 
 ## About
 
-`StructuralIdentifiability.jl` is a Julia package for assessing structural parameter identifiability of parametric ODE models, both local and global.
+`StructuralIdentifiability.jl` is a Julia package for assessing structural identifiability and observability of parametric ODE models, both local and global.
 This includes computation of identifiable functions of states and parameters. The package also offers functionality to assess local identifiability
 in discrete-time models.
 For an introduction to structural identifiability, we refer to [[3]](#review).
@@ -122,7 +122,7 @@ will return
 ```
 
 which are exactly the identifiable functions we have found before. Furthermore, by specifying `with_states = true`, one can compute the generating set for
-all identifiable functions of parameters and states (in other words, all observable functions):
+all identifiable functions of parameters and states (in other words, all observable functions, also known as the observation field):
 
 ```julia
 find_identifiable_functions(ode, with_states = true)
diff --git a/docs/src/assets/Project.toml b/docs/src/assets/Project.toml
index 30e21dee0..359116e04 100644
--- a/docs/src/assets/Project.toml
+++ b/docs/src/assets/Project.toml
@@ -1,11 +1,11 @@
 [deps]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 StructuralIdentifiability = "220ca800-aa68-49bb-acd8-6037fa93a544"
+ModelingToolkitBase = "7771a370-6774-4173-bd38-47e70ca0b839"
 
 [compat]
 BenchmarkTools = "1.3"
 Documenter = "0.27, 1"
-ModelingToolkit = "9"
+ModelingToolkitBase = "1.1"
 StructuralIdentifiability = "0.5"
diff --git a/docs/src/identifiability/identifiability.md b/docs/src/identifiability/identifiability.md
index cb172fd15..1c8752c73 100644
--- a/docs/src/identifiability/identifiability.md
+++ b/docs/src/identifiability/identifiability.md
@@ -17,3 +17,7 @@ assess_local_identifiability
 ```@docs
 find_identifiable_functions
 ```
+
+```@docs
+observation_field
+```
diff --git a/docs/src/index.md b/docs/src/index.md
index baddd7a3b..bbc9e2c61 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -28,13 +28,13 @@ Pkg.add("StructuralIdentifiability")
 
 ## Feature Summary
 
-`StructuralIdentifiability.jl` can assess local and global identifiability of ODE models. In addition to these straightforward identifiability queries on individual parameters, the package can distinguish between single- and multi-experiment identifiability.
+`StructuralIdentifiability.jl` can assess local and global identifiability of ODE models. In addition to these straightforward identifiability queries on individual parameters, the package can distinguish between single- and multi-experiment identifiability, and can compute identifiable functions of parameters and states.
 
 ## Feature List
 
   - Local identifiability checks
   - Global identifiability checks
-  - Assessment of identifiable functions of parameters and states
+  - Identifiable functions of parameters and states (a.k.a. observation field)
   - Model reparametrization (experimental)
 
 ## External tutorials
diff --git a/docs/src/tutorials/identifiable_functions.md b/docs/src/tutorials/identifiable_functions.md
index b297e6336..0598675bc 100644
--- a/docs/src/tutorials/identifiable_functions.md
+++ b/docs/src/tutorials/identifiable_functions.md
@@ -1,4 +1,4 @@
-# Globally Identifiable Functions
+# Globally Identifiable Functions (a.k.a. Observation Field)
 
 In addition to assessing identifiability of given functions of parameters or states, `StructuralIdentifiability.jl`
 provides the function `find_identifiable_functions` which finds a set of identifiable functions such that any other
@@ -30,6 +30,17 @@ of numbers `{p1, p3}` is uniquely determined since `p1 + p3` and `p1 * p3` are k
 Furthermore, we see that, for fixed input and output, `p2` and `p4` can take infinitely many values but
 knowing one of them, we would also be able to determine the other.
 
+!!! note "Observation field"
+    The field of identifiable functions involving both parameters and states is also called the **observation field**. It can be computed using
+    [`observation_field`](@ref):
+
+    ```julia
+    observation_field(LLW1987)
+    ```
+
+    This is equivalent to calling
+    `find_identifiable_functions(LLW1987, with_states = true)`.
+
 Moreover, we can find generators of all identifiable functions in parameters and states:
 
 ```@example funcs
diff --git a/src/StructuralIdentifiability.jl b/src/StructuralIdentifiability.jl
index 6aa9fedbf..67496f11e 100644
--- a/src/StructuralIdentifiability.jl
+++ b/src/StructuralIdentifiability.jl
@@ -95,6 +95,9 @@ export set_parameter_values, rename_variables
 # extra functionality
 export find_ioequations, find_identifiable_functions, default_cmp
 
+# observation field!!
+export observation_field
+
 # exporting to other formats
 export print_for_maple, print_for_DAISY, print_for_COMBOS, print_for_GenSSI
 
@@ -126,6 +129,7 @@ include("states.jl")
 include("output_saturation.jl")
 include("global_identifiability.jl")
 include("identifiable_functions.jl")
+include("observation_field.jl")
 include("parametrizations.jl")
 include("lincomp.jl")
 include("pb_representation.jl")
diff --git a/src/identifiable_functions.jl b/src/identifiable_functions.jl
index a186e0c64..b3d119e4b 100644
--- a/src/identifiable_functions.jl
+++ b/src/identifiable_functions.jl
@@ -31,21 +31,18 @@ This function takes the following optional arguments:
 ## Example
 
 ```jldoctest; setup = :(using Logging; Logging.disable_logging(Logging.Info);)
-using StructuralIdentifiability
+julia> using StructuralIdentifiability
   
-ode = @ODEmodel(
-    x0'(t) = -(a01 + a21) * x0(t) + a12 * x1(t),
-    x1'(t) = a21 * x0(t) - a12 * x1(t),
-    y(t) = x0(t)
-)
+julia> ode = @ODEmodel(
+           x0'(t) = -(a01 + a21) * x0(t) + a12 * x1(t),
+           x1'(t) = a21 * x0(t) - a12 * x1(t),
+           y(t) = x0(t)
+       );
 
-find_identifiable_functions(ode)
+julia> identifiable_functions = find_identifiable_functions(ode; loglevel = Logging.Error);
 
-# output
-
-2-element Vector{AbstractAlgebra.Generic.FracFieldElem{Nemo.QQMPolyRingElem}}:
- a01 + a12 + a21
- a01*a12
+julia> length(identifiable_functions)
+2
 ```
 
 """
diff --git a/src/observation_field.jl b/src/observation_field.jl
new file mode 100644
index 000000000..a4cda7923
--- /dev/null
+++ b/src/observation_field.jl
@@ -0,0 +1,49 @@
+"""
+    observation_field(ode::ODE; options...)
+
+Find simple generators of the observation field of the given ODE system.
+
+## Options
+
+This function takes the following optional arguments:
+- `known_ic`: a list of functions whose initial conditions are assumed to be known,
+  then the returned identifiable functions will be functions of parameters and
+  initial conditions, not states (this is an experimental functionality).
+- `prob_threshold`: A float in the range from 0 to 1, the probability of correctness. Default
+  is `0.99`.
+- `seed`: The rng seed. Default value is `42`.
+- `cmp`: A comparator for rational functions. For two rational functions `f1`
+    and `f2`, `cmp(f1, f2)` should be `true` if `f1` is simpler than `f2`.
+    By default, `StructuralIdentifiability.default_cmp(ode)` is used.
+- `loglevel` - the minimal level of log messages to display (`Logging.Info` by default)
+
+## Example
+
+```jldoctest; setup = :(using Logging; Logging.disable_logging(Logging.Info);)
+julia> using StructuralIdentifiability
+
+julia> ode = @ODEmodel(
+           x0'(t) = -(a01 + a21) * x0(t) + a12 * x1(t),
+           x1'(t) = a21 * x0(t) - a12 * x1(t),
+           y(t) = x0(t)
+       );
+
+julia> observation_field(ode; loglevel=Logging.Warn)
+4-element Vector{AbstractAlgebra.Generic.FracFieldElem{Nemo.QQMPolyRingElem}}:
+ a01 + a12 + a21
+ a01*a12
+ x0(t)
+ x0(t)*a12 + x1(t)*a12
+```
+
+"""
+function observation_field(
+        ode::ODE{T};
+        known_ic::Vector{<:ExtendedFraction{T}} = Vector{ExtendedFraction{T}}(),
+        prob_threshold::Float64 = 0.99,
+        seed = 42,
+        cmp = default_cmp(ode),
+        loglevel = Logging.Info,
+    ) where {T <: MPolyRingElem{QQFieldElem}}
+    return find_identifiable_functions(ode, with_states = true, known_ic = known_ic, prob_threshold = prob_threshold, seed = seed, loglevel = loglevel, cmp = cmp)
+end
diff --git a/test/bodies/observation_field.jl b/test/bodies/observation_field.jl
new file mode 100644
index 000000000..0e86a127f
--- /dev/null
+++ b/test/bodies/observation_field.jl
@@ -0,0 +1,11 @@
+include(joinpath(@__DIR__, "..", "shared", "test_setup.jl"))
+
+@testset "Identifiable functions of parameters" begin
+    ode = @ODEmodel(
+        x1'(t) = -(k21 + V_M / (K_M + x1(t))) * x1(t) + k12 * x2(t) + b1 * u(t),
+        x2'(t) = k21 * x1(t) - (k02 + k12) * x2(t),
+        y(t) = c * x1(t)
+    )
+
+    @test observation_field(ode) == find_identifiable_functions(ode, with_states = true)
+end