diff --git a/.github/workflows/PerformanceCheck.yml b/.github/workflows/PerformanceCheck.yml
new file mode 100644
index 000000000..2ed2c5eef
--- /dev/null
+++ b/.github/workflows/PerformanceCheck.yml
@@ -0,0 +1,29 @@
+name: Performance Check
+on:
+  push:
+    branches:
+      - 'master'
+    tags: '*'
+  pull_request:
+jobs:
+  build-and-benchmark:
+    runs-on: ${{ matrix.os }}
+    strategy:
+      matrix:
+        julia-version: ['1']
+        julia-arch: [x64]
+        os: [ubuntu-latest]
+    steps:
+    - name: Checkout
+      uses: actions/checkout@v2
+    - name: Install Julia
+      uses: julia-actions/setup-julia@latest
+      with:
+        version: ${{ matrix.julia-version }}
+    - name: Build project
+      uses: julia-actions/julia-buildpkg@latest
+    - name: Benchmark
+      env:
+        GITHUB_TOKEN: ${{ secrets.BENCHMARK_KEY }}
+        PR_NUMBER: ${{ github.event.number }}
+      run: julia --threads=4 --project=./benchmarking/CI-scripts ./benchmarking/CI-scripts/runtests.jl 2
diff --git a/benchmarking/CI-scripts/Project.toml b/benchmarking/CI-scripts/Project.toml
new file mode 100644
index 000000000..2f6685d27
--- /dev/null
+++ b/benchmarking/CI-scripts/Project.toml
@@ -0,0 +1,8 @@
+[deps]
+ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
+GitHub = "bc5e4493-9b4d-5f90-b8aa-2b2bcaad7a26"
+GitHubActions = "6b79fd1a-b13a-48ab-b6b0-aaee1fee41df"
+Nemo = "2edaba10-b0f1-5616-af89-8c11ac63239a"
+PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
+Primes = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae"
+TestSetExtensions = "98d24dd4-01ad-11ea-1b02-c9a08f80db04"
diff --git a/benchmarking/CI-scripts/benchmarks.jl b/benchmarking/CI-scripts/benchmarks.jl
new file mode 100644
index 000000000..a5b2afba6
--- /dev/null
+++ b/benchmarking/CI-scripts/benchmarks.jl
@@ -0,0 +1,450 @@
+benchmarks = Dict(
+    :LV_simple => Dict(
+        :name => "Modified LV for testing",
+        :ode => @ODEmodel(
+            x1'(t) = (a + b) * x1(t) - c * x1(t) * x2(t),
+            x2'(t) = -a * b * x2(t) + d * x1(t) * x2(t),
+            y1(t) = x1(t)
+        ),
+    ),
+    :SIWR => Dict(
+        :name => "SIWR original",
+        :ode => @ODEmodel(
+            S'(t) = mu - bi * S(t) * I(t) - bw * S(t) * W(t) - mu * S(t) + a * R(t),
+            I'(t) = bw * S(t) * W(t) + bi * S(t) * I(t) - (gam + mu) * I(t),
+            W'(t) = xi * (I(t) - W(t)),
+            R'(t) = gam * I(t) - (mu + a) * R(t),
+            y(t) = k * I(t)
+        ),
+        :skip => false,
+    ),
+    :SIWR_simplified => Dict(
+        :name => "SIWR with extra output",
+        :ode => @ODEmodel(
+            S'(t) = mu - bi * S(t) * I(t) - bw * S(t) * W(t) - mu * S(t) + a * R(t),
+            I'(t) = bw * S(t) * W(t) + bi * S(t) * I(t) - (gam + mu) * I(t),
+            W'(t) = xi * (I(t) - W(t)),
+            R'(t) = gam * I(t) - (mu + a) * R(t),
+            y(t) = k * I(t),
+            y2(t) = S(t) + I(t) + R(t)
+        ),
+        :skip => false,
+    ),
+    :Pharm => Dict(
+        :name => "Pharm",
+        :ode => @ODEmodel(
+            x0'(t) =
+                a1 * (x1(t) - x0(t)) -
+                (ka * n * x0(t)) / (kc * ka + kc * x2(t) + ka * x0(t)),
+            x1'(t) = a2 * (x0(t) - x1(t)),
+            x2'(t) =
+                b1 * (x3(t) - x2(t)) -
+                (kc * n * x2(t)) / (kc * ka + kc * x2(t) + ka * x0(t)),
+            x3'(t) = b2 * (x2(t) - x3(t)),
+            y1(t) = x0(t)
+        ),
+        :skip => false,
+    ),
+    :SEAIJRC => Dict(
+        :name => "SEAIJRC Covid model",
+        :ode => @ODEmodel(
+            S'(t) = -b * S(t) * (I(t) + J(t) + q * A(t)) * Ninv,
+            E'(t) = b * S(t) * (I(t) + J(t) + q * A(t)) * Ninv - k * E(t),
+            A'(t) = k * (1 - r) * E(t) - g1 * A(t),
+            I'(t) = k * r * E(t) - (alpha + g1) * I(t),
+            J'(t) = alpha * I(t) - g2 * J(t),
+            C'(t) = alpha * I(t),
+            y(t) = C(t),
+            y2(t) = Ninv
+        ),
+        :skip => false,
+    ),
+    :MAPK5 => Dict(
+        :name => "MAPK model (5 outputs)",
+        :ode => @ODEmodel(
+            KS00'(t) =
+                -a00 * K(t) * S00(t) +
+                b00 * KS00(t) +
+                gamma0100 * FS01(t) +
+                gamma1000 * FS10(t) +
+                gamma1100 * FS11(t),
+            KS01'(t) =
+                -a01 * K(t) * S01(t) + b01 * KS01(t) + c0001 * KS00(t) -
+                alpha01 * F(t) * S01(t) +
+                beta01 * FS01(t) +
+                gamma1101 * FS11(t),
+            KS10'(t) =
+                -a10 * K(t) * S10(t) + b10 * KS10(t) + c0010 * KS00(t) -
+                alpha10 * F(t) * S10(t) +
+                beta10 * FS10(t) +
+                gamma1110 * FS11(t),
+            FS01'(t) =
+                -alpha11 * F(t) * S11(t) +
+                beta11 * FS11(t) +
+                c0111 * KS01(t) +
+                c1011 * KS10(t) +
+                c0011 * KS00(t),
+            FS10'(t) = a00 * K(t) * S00(t) - (b00 + c0001 + c0010 + c0011) * KS00(t),
+            FS11'(t) = a01 * K(t) * S01(t) - (b01 + c0111) * KS01(t),
+            K'(t) = a10 * K(t) * S10(t) - (b10 + c1011) * KS10(t),
+            F'(t) = alpha01 * F(t) * S01(t) - (beta01 + gamma0100) * FS01(t),
+            S00'(t) = alpha10 * F(t) * S10(t) - (beta10 + gamma1000) * FS10(t),
+            S01'(t) =
+                alpha11 * F(t) * S11(t) -
+                (beta11 + gamma1101 + gamma1110 + gamma1100) * FS11(t),
+            S10'(t) =
+                -a00 * K(t) * S00(t) + (b00 + c0001 + c0010 + c0011) * KS00(t) -
+                a01 * K(t) * S01(t) + (b01 + c0111) * KS01(t) -
+                a10 * K(t) * S10(t) + (b10 + c1011) * KS10(t),
+            S11'(t) =
+                -alpha01 * F(t) * S01(t) + (beta01 + gamma0100) * FS01(t) -
+                alpha10 * F(t) * S10(t) + (beta10 + gamma1000) * FS10(t) -
+                alpha11 * F(t) * S11(t) +
+                (beta11 + gamma1101 + gamma1110 + gamma1100) * FS11(t),
+            y1(t) = F(t),
+            y2(t) = S00(t),
+            y3(t) = S01(t),
+            y4(t) = S10(t),
+            y5(t) = S11(t)
+        ),
+    ),
+    :MAPK5bis => Dict(
+        :name => "MAPK model (5 outputs bis)",
+        :ode => @ODEmodel(
+            KS00'(t) =
+                -a00 * K(t) * S00(t) +
+                b00 * KS00(t) +
+                gamma0100 * FS01(t) +
+                gamma1000 * FS10(t) +
+                gamma1100 * FS11(t),
+            KS01'(t) =
+                -a01 * K(t) * S01(t) + b01 * KS01(t) + c0001 * KS00(t) -
+                alpha01 * F(t) * S01(t) +
+                beta01 * FS01(t) +
+                gamma1101 * FS11(t),
+            KS10'(t) =
+                -a10 * K(t) * S10(t) + b10 * KS10(t) + c0010 * KS00(t) -
+                alpha10 * F(t) * S10(t) +
+                beta10 * FS10(t) +
+                gamma1110 * FS11(t),
+            FS01'(t) =
+                -alpha11 * F(t) * S11(t) +
+                beta11 * FS11(t) +
+                c0111 * KS01(t) +
+                c1011 * KS10(t) +
+                c0011 * KS00(t),
+            FS10'(t) = a00 * K(t) * S00(t) - (b00 + c0001 + c0010 + c0011) * KS00(t),
+            FS11'(t) = a01 * K(t) * S01(t) - (b01 + c0111) * KS01(t),
+            K'(t) = a10 * K(t) * S10(t) - (b10 + c1011) * KS10(t),
+            F'(t) = alpha01 * F(t) * S01(t) - (beta01 + gamma0100) * FS01(t),
+            S00'(t) = alpha10 * F(t) * S10(t) - (beta10 + gamma1000) * FS10(t),
+            S01'(t) =
+                alpha11 * F(t) * S11(t) -
+                (beta11 + gamma1101 + gamma1110 + gamma1100) * FS11(t),
+            S10'(t) =
+                -a00 * K(t) * S00(t) + (b00 + c0001 + c0010 + c0011) * KS00(t) -
+                a01 * K(t) * S01(t) + (b01 + c0111) * KS01(t) -
+                a10 * K(t) * S10(t) + (b10 + c1011) * KS10(t),
+            S11'(t) =
+                -alpha01 * F(t) * S01(t) + (beta01 + gamma0100) * FS01(t) -
+                alpha10 * F(t) * S10(t) + (beta10 + gamma1000) * FS10(t) -
+                alpha11 * F(t) * S11(t) +
+                (beta11 + gamma1101 + gamma1110 + gamma1100) * FS11(t),
+            y0(t) = K(t),
+            y1(t) = F(t),
+            y2(t) = S00(t),
+            y3(t) = S01(t) + S10(t),
+            y4(t) = S11(t)
+        ),
+        :skip => false,
+    ),
+    :MAPK6 => Dict(
+        :name => "MAPK model (6 outputs)",
+        :ode => @ODEmodel(
+            KS00'(t) =
+                -a00 * K(t) * S00(t) +
+                b00 * KS00(t) +
+                gamma0100 * FS01(t) +
+                gamma1000 * FS10(t) +
+                gamma1100 * FS11(t),
+            KS01'(t) =
+                -a01 * K(t) * S01(t) + b01 * KS01(t) + c0001 * KS00(t) -
+                alpha01 * F(t) * S01(t) +
+                beta01 * FS01(t) +
+                gamma1101 * FS11(t),
+            KS10'(t) =
+                -a10 * K(t) * S10(t) + b10 * KS10(t) + c0010 * KS00(t) -
+                alpha10 * F(t) * S10(t) +
+                beta10 * FS10(t) +
+                gamma1110 * FS11(t),
+            FS01'(t) =
+                -alpha11 * F(t) * S11(t) +
+                beta11 * FS11(t) +
+                c0111 * KS01(t) +
+                c1011 * KS10(t) +
+                c0011 * KS00(t),
+            FS10'(t) = a00 * K(t) * S00(t) - (b00 + c0001 + c0010 + c0011) * KS00(t),
+            FS11'(t) = a01 * K(t) * S01(t) - (b01 + c0111) * KS01(t),
+            K'(t) = a10 * K(t) * S10(t) - (b10 + c1011) * KS10(t),
+            F'(t) = alpha01 * F(t) * S01(t) - (beta01 + gamma0100) * FS01(t),
+            S00'(t) = alpha10 * F(t) * S10(t) - (beta10 + gamma1000) * FS10(t),
+            S01'(t) =
+                alpha11 * F(t) * S11(t) -
+                (beta11 + gamma1101 + gamma1110 + gamma1100) * FS11(t),
+            S10'(t) =
+                -a00 * K(t) * S00(t) + (b00 + c0001 + c0010 + c0011) * KS00(t) -
+                a01 * K(t) * S01(t) + (b01 + c0111) * KS01(t) -
+                a10 * K(t) * S10(t) + (b10 + c1011) * KS10(t),
+            S11'(t) =
+                -alpha01 * F(t) * S01(t) + (beta01 + gamma0100) * FS01(t) -
+                alpha10 * F(t) * S10(t) + (beta10 + gamma1000) * FS10(t) -
+                alpha11 * F(t) * S11(t) +
+                (beta11 + gamma1101 + gamma1110 + gamma1100) * FS11(t),
+            y0(t) = K(t),
+            y1(t) = F(t),
+            y2(t) = S00(t),
+            y3(t) = S01(t),
+            y4(t) = S10(t),
+            y5(t) = S11(t)
+        ),
+    ),
+    :Goodwin => Dict(
+        :name => "Goodwin oscillator",
+        :ode => @ODEmodel(
+            x1'(t) = -b * x1(t) + 1 / (c + x4(t)),
+            x2'(t) = alpha * x1(t) - beta * x2(t),
+            x3'(t) = gama * x2(t) - delta * x3(t),
+            x4'(t) = sigma * x4(t) * (gama * x2(t) - delta * x3(t)) / x3(t),
+            y(t) = x1(t)
+        ),
+    ),
+    :HIV => Dict(
+        :name => "HIV",
+        :ode => @ODEmodel(
+            x'(t) = lm - d * x(t) - beta * x(t) * v(t),
+            y'(t) = beta * x(t) * v(t) - a * y(t),
+            v'(t) = k * y(t) - u * v(t),
+            w'(t) = c * x(t) * y(t) * w(t) - c * q * y(t) * w(t) - b * w(t),
+            z'(t) = c * q * y(t) * w(t) - h * z(t),
+            y1(t) = w(t),
+            y2(t) = z(t)
+        ),
+    ),
+    :SIRSforced => Dict(
+        :name => "SIRS forced",
+        :ode => @ODEmodel(
+            s'(t) = mu - mu * s(t) - b0 * (1 + b1 * x1(t)) * i(t) * s(t) + g * r(t),
+            i'(t) = b0 * (1 + b1 * x1(t)) * i(t) * s(t) - (nu + mu) * i(t),
+            r'(t) = nu * i(t) - (mu + g) * r(t),
+            x1'(t) = -M * x2(t),
+            x2'(t) = M * x1(t),
+            y1(t) = i(t),
+            y2(t) = r(t)
+        ),
+    ),
+    :Akt => Dict(
+        :name => "Akt pathway",
+        :ode => @ODEmodel(
+            EGFR'(t) =
+                EGFR_turnover * pro_EGFR(t) + EGF_EGFR(t) * reaction_1_k2 -
+                EGFR(t) * EGFR_turnover - EGF_EGFR(t) * reaction_1_k1,
+            pEGFR'(t) =
+                EGF_EGFR(t) * reaction_9_k1 - pEGFR(t) * reaction_4_k1 +
+                pEGFR_Akt(t) * reaction_2_k2 +
+                pEGFR_Akt(t) * reaction_3_k1 -
+                Akt(t) * pEGFR(t) * reaction_2_k1,
+            pEGFR_Akt'(t) =
+                Akt(t) * pEGFR(t) * reaction_2_k1 - pEGFR_Akt(t) * reaction_3_k1 - pEGFR_Akt(t) * reaction_2_k2,
+            Akt'(t) =
+                pAkt(t) * reaction_7_k1 + pEGFR_Akt(t) * reaction_2_k2 -
+                Akt(t) * pEGFR(t) * reaction_2_k1,
+            pAkt'(t) =
+                pAkt_S6(t) * reaction_5_k2 - pAkt(t) * reaction_7_k1 +
+                pAkt_S6(t) * reaction_6_k1 +
+                pEGFR_Akt(t) * reaction_3_k1 - S6(t) * pAkt(t) * reaction_5_k1,
+            S6'(t) =
+                pAkt_S6(t) * reaction_5_k2 + pS6(t) * reaction_8_k1 -
+                S6(t) * pAkt(t) * reaction_5_k1,
+            pAkt_S6'(t) =
+                S6(t) * pAkt(t) * reaction_5_k1 - pAkt_S6(t) * reaction_6_k1 -
+                pAkt_S6(t) * reaction_5_k2,
+            pS6'(t) = pAkt_S6(t) * reaction_6_k1 - pS6(t) * reaction_8_k1,
+            EGF_EGFR'(t) =
+                EGF_EGFR(t) * reaction_1_k1 - EGF_EGFR(t) * reaction_9_k1 -
+                EGF_EGFR(t) * reaction_1_k2,
+            y1(t) = a1 * (pEGFR(t) + pEGFR_Akt(t)),
+            y2(t) = a2 * (pAkt(t) + pAkt_S6(t)),
+            y3(t) = a3 * pS6(t)
+        ),
+    ),
+    :CD8_Tcell => Dict(
+        :name => "CD8 T cell differentiation",
+        :ode => @ODEmodel(
+            N'(t) = -N(t) * mu_N - N(t) * P(t) * delta_NE,
+            E'(t) =
+                N(t) * P(t) * delta_NE - E(t)^2 * mu_EE - E(t) * delta_EL +
+                E(t) * P(t) * rho_E,
+            S'(t) =
+                S(t) * delta_EL - S(t) * delta_LM - S(t)^2 * mu_LL - E(t) * S(t) * mu_LE,
+            M'(t) = S(t) * delta_LM - mu_M * M(t),
+            P'(t) =
+                P(t)^2 * rho_P - P(t) * mu_P - E(t) * P(t) * mu_PE - S(t) * P(t) * mu_PL,
+            y1(t) = N(t),
+            y2(t) = E(t) + S(t),
+            y3(t) = M(t)
+        ),
+    ),
+    :CRN => Dict(
+        :name => "Chemical reaction network",
+        :ode => @ODEmodel(
+            x1'(t) = -k1 * x1(t) * x2(t) + k2 * x4(t) + k4 * x6(t),
+            x2'(t) = -k1 * x1(t) * x2(t) + k2 * x4(t) + k3 * x4(t),
+            x3'(t) = k3 * x4(t) + k5 * x6(t) - k6 * x3(t) * x5(t),
+            x4'(t) = k1 * x1(t) * x2(t) - k2 * x4(t) - k3 * x4(t),
+            x5'(t) = k4 * x6(t) + k5 * x6(t) - k6 * x3(t) * x5(t),
+            x6'(t) = -k4 * x6(t) - k5 * x6(t) + k6 * x3(t) * x5(t),
+            y1(t) = x3(t),
+            y2(t) = x2(t)
+        ),
+    ),
+    :SLIQR => Dict(
+        :name => "SLIQR",
+        :ode => @ODEmodel(
+            S'(t) = -b * In(t) * S(t) * Ninv - u(t) * S(t) * Ninv,
+            L'(t) = b * In(t) * S(t) * Ninv - a * L(t),
+            In'(t) = a * L(t) - g * In(t) + s * Q(t),
+            Q'(t) = (1 - e) * g * In(t) - s * Q(t),
+            y(t) = In(t) * Ninv
+        ),
+    ),
+    :Sntg => Dict(
+        :name => "Sntg",
+        :ode => @ODEmodel(
+            S'(t) = r * S(t) - (e + a * W(t)) * S(t) - d * W(t) * S(t) + g * R(t),
+            R'(t) = rR * R(t) + (e + a * W(t)) * S(t) - dr * W(t) * R(t) - g * R(t),
+            W'(t) = Dd * (T - W(t)),
+            y1(t) = S(t) + R(t),
+            y2(t) = T
+        ),
+    ),
+    :QY => Dict(
+        :name => "QY",
+        :ode => @ODEmodel(
+            P0'(t) = P1(t),
+            P1'(t) = P2(t),
+            P2'(t) = P3(t),
+            P3'(t) = P4(t),
+            P4'(t) =
+                -(
+                    Ks * M * siga1 * siga2 * P1(t) +
+                    (
+                        Ks * M * siga1 +
+                        Ks * M * siga2 +
+                        Ks * siga1 * siga2 +
+                        siga1 * siga2 * M
+                    ) * P2(t) +
+                    (
+                        Ks * M +
+                        Ks * siga1 +
+                        Ks * siga2 +
+                        M * siga1 +
+                        M * siga2 +
+                        siga1 * siga2
+                    ) * P3(t) +
+                    (Ks + M + siga1 + siga2) * P4(t)
+                ) -
+                (
+                    Mar * P5(t) +
+                    beta +
+                    beta_SA / (siga2 * M) * (
+                        P3(t) +
+                        P2(t) * (Ks + M + Mar) +
+                        P1(t) * (Ks * M + Ks * Mar + M * Mar) +
+                        P0(t) * Ks * M * Mar
+                    ) +
+                    beta_SI / M * (P2(t) + P1(t) * (Ks + Mar) + P0(t) * Ks * Mar) +
+                    beta_SA * phi / ((1 - phi) * siga2 * M) * (
+                        P3(t) +
+                        P2(t) * (Ks + M + siga2) +
+                        P1(t) * (Ks * M + Ks * siga2 + M * siga2) +
+                        P0(t) * Ks * M * siga2
+                    )
+                ) * (
+                    alpa +
+                    Ks * M * siga1 * siga2 * P0(t) +
+                    (
+                        Ks * M * siga1 +
+                        Ks * M * siga2 +
+                        Ks * siga1 * siga2 +
+                        siga1 * siga2 * M
+                    ) * P1(t) +
+                    (
+                        Ks * M +
+                        Ks * siga1 +
+                        Ks * siga2 +
+                        M * siga1 +
+                        M * siga2 +
+                        siga1 * siga2
+                    ) * P2(t) +
+                    (Ks + M + siga1 + siga2) * P3(t) +
+                    P4(t)
+                ),
+            P5'(t) =
+                -Mar * P5(t) - (
+                    beta +
+                    beta_SA / (siga2 * M) * (
+                        P3(t) +
+                        P2(t) * (Ks + M + Mar) +
+                        P1(t) * (Ks * M + Ks * Mar + M * Mar) +
+                        P0(t) * Ks * M * Mar
+                    ) +
+                    beta_SI / M * (P2(t) + P1(t) * (Ks + Mar) + P0(t) * Ks * Mar) +
+                    beta_SA * phi / ((1 - phi) * siga2 * M) * (
+                        P3(t) +
+                        P2(t) * (Ks + M + siga2) +
+                        P1(t) * (Ks * M + Ks * siga2 + M * siga2) +
+                        P0(t) * Ks * M * siga2
+                    )
+                ),
+            y(t) = P0(t)
+        ),
+    ),
+    :LLW => Dict(
+        # https://github.com/Xabo-RB/Local-Global-Models/blob/main/Models/General/LLW1987_io.jl
+        :name => "LLW1987_io",
+        :ode => @ODEmodel(
+            x1'(t) = -p1 * x1(t) + p2 * u(t),
+            x2'(t) = -p3 * x2(t) + p4 * u(t),
+            x3'(t) = -(p1 + p3) * x3(t) + (p4 * x1(t) + p2 * x2(t)) * u(t),
+            y1(t) = x3(t)
+        ),
+    ),
+    :Bilirubin => Dict(
+        # https://github.com/Xabo-RB/Local-Global-Models/blob/main/Models/Physiology/Bilirubin2_io.jl
+        :name => "Bilirubin2_io",
+        :ode => @ODEmodel(
+            x1'(t) =
+                -(k21 + k31 + k41 + k01) * x1(t) +
+                k12 * x2(t) +
+                k13 * x3(t) +
+                k14 * x4(t) +
+                u(t),
+            x2'(t) = k21 * x1(t) - k12 * x2(t),
+            x3'(t) = k31 * x1(t) - k13 * x3(t),
+            x4'(t) = k41 * x1(t) - k14 * x4(t),
+            y1(t) = x1(t)
+        ),
+    ),
+    # https://github.com/Xabo-RB/Local-Global-Models/blob/main/Models/Epidemiology/SEUIR.jl
+    :SEUIR => Dict(
+        :name => "SEUIR",
+        :ode => @ODEmodel(
+            S'(t) = -beta * (U(t) + I(t)) * (S(t) / N),
+            E'(t) = beta * (U(t) + I(t)) * (S(t) / N) - E(t) * z,
+            U'(t) = (z - w) * E(t) - U(t) * d,
+            I'(t) = w * E(t) - I(t) * d,
+            R'(t) = (U(t) + I(t)) * d,
+            y1(t) = I(t)
+        ),
+    ),
+)
diff --git a/benchmarking/CI-scripts/common.jl b/benchmarking/CI-scripts/common.jl
new file mode 100644
index 000000000..b9d91d65e
--- /dev/null
+++ b/benchmarking/CI-scripts/common.jl
@@ -0,0 +1,12 @@
+
+function dump_results(file, key)
+    open(file, "w") do out
+        println(out, key)
+        for problem in suite
+            problem_name = problem.problem_name
+            result = problem.result
+            type = problem.type
+            println(out, "$problem_name:$(String(type)):$(join(map(string, result), ","))")
+        end
+    end
+end
diff --git a/benchmarking/CI-scripts/run-on-master/Project.toml b/benchmarking/CI-scripts/run-on-master/Project.toml
new file mode 100644
index 000000000..2f6685d27
--- /dev/null
+++ b/benchmarking/CI-scripts/run-on-master/Project.toml
@@ -0,0 +1,8 @@
+[deps]
+ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
+GitHub = "bc5e4493-9b4d-5f90-b8aa-2b2bcaad7a26"
+GitHubActions = "6b79fd1a-b13a-48ab-b6b0-aaee1fee41df"
+Nemo = "2edaba10-b0f1-5616-af89-8c11ac63239a"
+PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
+Primes = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae"
+TestSetExtensions = "98d24dd4-01ad-11ea-1b02-c9a08f80db04"
diff --git a/benchmarking/CI-scripts/run-on-master/run_benchmarks.jl b/benchmarking/CI-scripts/run-on-master/run_benchmarks.jl
new file mode 100644
index 000000000..e3f97cab3
--- /dev/null
+++ b/benchmarking/CI-scripts/run-on-master/run_benchmarks.jl
@@ -0,0 +1,13 @@
+using Pkg
+Pkg.activate(@__DIR__)
+Pkg.instantiate()
+Pkg.status()
+
+Pkg.add(url = "https://github.com/SciML/StructuralIdentifiability.jl")
+
+id = ARGS[1]
+
+include("../common.jl")
+include("../run_benchmarks.jl")
+
+dump_results((@__DIR__) * "/results_$id", "master")
diff --git a/benchmarking/CI-scripts/run-on-nightly/Project.toml b/benchmarking/CI-scripts/run-on-nightly/Project.toml
new file mode 100644
index 000000000..2f6685d27
--- /dev/null
+++ b/benchmarking/CI-scripts/run-on-nightly/Project.toml
@@ -0,0 +1,8 @@
+[deps]
+ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
+GitHub = "bc5e4493-9b4d-5f90-b8aa-2b2bcaad7a26"
+GitHubActions = "6b79fd1a-b13a-48ab-b6b0-aaee1fee41df"
+Nemo = "2edaba10-b0f1-5616-af89-8c11ac63239a"
+PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
+Primes = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae"
+TestSetExtensions = "98d24dd4-01ad-11ea-1b02-c9a08f80db04"
diff --git a/benchmarking/CI-scripts/run-on-nightly/run_benchmarks.jl b/benchmarking/CI-scripts/run-on-nightly/run_benchmarks.jl
new file mode 100644
index 000000000..69032b41d
--- /dev/null
+++ b/benchmarking/CI-scripts/run-on-nightly/run_benchmarks.jl
@@ -0,0 +1,12 @@
+import Pkg
+Pkg.activate(@__DIR__)
+Pkg.develop(path = (@__DIR__) * "/../../../")
+Pkg.instantiate()
+Pkg.status()
+
+id = ARGS[1]
+
+include("../common.jl")
+include("../run_benchmarks.jl")
+
+dump_results((@__DIR__) * "/results_$id", "nightly")
diff --git a/benchmarking/CI-scripts/run_benchmarks.jl b/benchmarking/CI-scripts/run_benchmarks.jl
new file mode 100644
index 000000000..5001f3c3e
--- /dev/null
+++ b/benchmarking/CI-scripts/run_benchmarks.jl
@@ -0,0 +1,127 @@
+# Adapter from https://github.com/sumiya11/Groebner.jl/blob/master/benchmark/CI-scripts/runtests.jl
+# GPL license
+
+# This script is not to be run directly. See runtests.jl.
+
+stopwatch = time_ns()
+
+# TTFX
+t1 = @timed using StructuralIdentifiability
+
+include("benchmarks.jl")
+
+suite = []
+
+push!(
+    suite,
+    (problem_name = "using StructuralIdentifiability", type = :time, result = [t1.time]),
+)
+
+# Allocations
+# Coming soon
+
+# Runtime
+import Primes
+import Nemo
+
+function nemo_make_prime_finite_field(p)
+    if p < typemax(UInt)
+        Nemo.fpField(convert(UInt, p), false)
+    else
+        Nemo.FpField(Nemo.ZZRingElem(p), false)
+    end
+end
+
+function perform_operation(
+    ode,
+    fun = StructuralIdentifiability.assess_identifiability,
+    trials = 3;
+    kws...,
+)
+    @assert fun in [
+        StructuralIdentifiability.assess_identifiability,
+        StructuralIdentifiability.assess_local_identifiability,
+        StructuralIdentifiability.find_identifiable_functions,
+        StructuralIdentifiability.reparametrize_global,
+    ]
+    times = []
+    for _ in 1:trials
+        GC.gc()
+        time = @elapsed fun(ode; kws...)
+        push!(times, time)
+    end
+    times
+end
+
+fun_name = "assess_identifiability"
+perform_operation(benchmarks[:LV_simple][:ode])
+for m in [:SIWR, :LV_simple, :Pharm]
+    push!(
+        suite,
+        (
+            problem_name = benchmarks[m][:name] * " " * fun_name,
+            type = :time,
+            result = perform_operation(benchmarks[m][:ode]),
+        ),
+    )
+end
+
+fun_name = "assess_local_identifiability"
+perform_operation(
+    benchmarks[:LV_simple][:ode],
+    StructuralIdentifiability.assess_local_identifiability,
+)
+for m in [:SIWR, :Pharm, :MAPK5, :MAPK5bis, :Goodwin]
+    push!(
+        suite,
+        (
+            problem_name = benchmarks[m][:name] * " " * fun_name,
+            type = :time,
+            result = perform_operation(
+                benchmarks[m][:ode],
+                StructuralIdentifiability.assess_local_identifiability,
+            ),
+        ),
+    )
+end
+
+fun_name = "find_identifiable_functions"
+perform_operation(
+    benchmarks[:LV_simple][:ode],
+    StructuralIdentifiability.find_identifiable_functions,
+)
+for m in [:Goodwin, :SEAIJRC, :Sntg, :LLW, :Bilirubin]
+    push!(
+        suite,
+        (
+            problem_name = benchmarks[m][:name] * " " * fun_name,
+            type = :time,
+            result = perform_operation(
+                benchmarks[m][:ode],
+                StructuralIdentifiability.find_identifiable_functions,
+            ),
+        ),
+    )
+end
+
+fun_name = "reparametrize_global"
+perform_operation(
+    benchmarks[:LV_simple][:ode],
+    StructuralIdentifiability.reparametrize_global,
+)
+for m in [:Goodwin, :SEAIJRC, :SEUIR]
+    push!(
+        suite,
+        (
+            problem_name = benchmarks[m][:name] * " " * fun_name,
+            type = :time,
+            result = perform_operation(
+                benchmarks[m][:ode],
+                StructuralIdentifiability.reparametrize_global,
+            ),
+        ),
+    )
+end
+
+stopwatch = time_ns() - stopwatch
+push!(suite, (problem_name = "total", type = :time, result = [stopwatch / 1e9]))
diff --git a/benchmarking/CI-scripts/runtests.jl b/benchmarking/CI-scripts/runtests.jl
new file mode 100644
index 000000000..a93a5d634
--- /dev/null
+++ b/benchmarking/CI-scripts/runtests.jl
@@ -0,0 +1,201 @@
+# Adapter from https://github.com/sumiya11/Groebner.jl/blob/master/benchmark/CI-scripts/runtests.jl
+# GPL license
+
+# Test for performance regressions.
+using Pkg
+Pkg.activate(@__DIR__)
+Pkg.instantiate()
+
+using ArgParse, GitHubActions, GitHub, Random, Logging
+using Test, TestSetExtensions, InteractiveUtils, PrettyTables
+using Base.Threads, Statistics
+
+const MAX_DEVIATION = 0.2
+const IGNORE_SMALL = 1e-3
+const SAMPLES = length(ARGS) > 0 ? parse(Int, ARGS[1]) : 1
+
+const dir_master = (@__DIR__) * "/run-on-master"
+const dir_nightly = (@__DIR__) * "/run-on-nightly"
+
+function runbench()
+    @info "Start benchmarking.."
+    @info "Using $(nthreads()) threads"
+    @info "Using $SAMPLES samples"
+
+    for i in 1:SAMPLES
+        # Run benchmarks on master
+        @info "Benchmarking StructuralIdentifiability.jl, master, running $dir_master"
+        @time run(
+            `$(Base.julia_cmd()) --startup-file=no --threads=$(nthreads()) --project=$dir_master $dir_master/run_benchmarks.jl $i`,
+            wait = true,
+        )
+
+        # Run benchmarks on nightly
+        @info "Benchmarking StructuralIdentifiability.jl, nightly, running $dir_nightly"
+        @time run(
+            `$(Base.julia_cmd()) --startup-file=no --threads=$(nthreads()) --project=$dir_nightly $dir_nightly/run_benchmarks.jl $i`,
+            wait = true,
+        )
+    end
+end
+
+# Adapted from https://github.com/MakieOrg/Makie.jl/blob/v0.21.0/metrics/ttfp/run-benchmark.jl.
+# License is MIT.
+function best_unit(m)
+    if m < 1e3
+        return 1, "ns"
+    elseif m < 1e6
+        return 1e3, "μs"
+    elseif m < 1e9
+        return 1e6, "ms"
+    else
+        return 1e9, "s"
+    end
+end
+
+function load_data()
+    results = []
+    for i in 1:SAMPLES
+        results_master_i = nothing
+        results_nightly_i = nothing
+        try
+            results_master_file = open(dir_master * "/results_$i", "r")
+            @info "Reading $results_master_file"
+            results_master_i = readlines(results_master_file)
+            close(results_master_file)
+        catch e
+            @warn "Error when reading the file with results"
+        end
+        try
+            results_nightly_file = open(dir_nightly * "/results_$i", "r")
+            @info "Reading $results_nightly_file"
+            results_nightly_i = readlines(results_nightly_file)
+            close(results_nightly_file)
+        catch e
+            @warn "Error when reading the file with results, sample $i"
+        end
+        @assert !(results_master_i === nothing) && !(results_nightly_i === nothing)
+        @assert length(results_master_i) == length(results_nightly_i)
+        @assert !isempty(results_master_i)
+        @assert results_master_i[1] == "master" && results_nightly_i[1] == "nightly"
+        results_master_i = results_master_i[2:end]
+        results_nightly_i = results_nightly_i[2:end]
+        push!(results, (master = results_master_i, nightly = results_nightly_i))
+    end
+    results
+end
+
+function clean_data(results)
+    nrecords = length(results[1][1])
+    results_problems = Vector{Any}(undef, nrecords)
+    results_types = Vector{Any}(undef, nrecords)
+    results_master = [[] for _ in 1:nrecords]
+    results_nightly = [[] for _ in 1:nrecords]
+    for i in 1:SAMPLES
+        for j in 1:nrecords
+            master = results[i].master[j]
+            nightly = results[i].nightly[j]
+            problem_name_master, type, times_master = split(master, ":")
+            problem_name_nightly, type, times_nightly = split(nightly, ":")
+            @assert problem_name_master == problem_name_nightly
+            times_master = map(
+                x -> parse(Float64, String(strip(x, ['[', ']', ' ', '\t']))),
+                split(times_master, ","),
+            )
+            times_nightly = map(
+                x -> parse(Float64, String(strip(x, ['[', ']', ' ', '\t']))),
+                split(times_nightly, ","),
+            )
+            append!(results_master[j], times_master)
+            append!(results_nightly[j], times_nightly)
+            results_problems[j] = join(split(problem_name_master, ","), " ")
+            results_types[j] = type
+        end
+    end
+    results_problems, results_types, results_master, results_nightly
+end
+
+# Compare results
+function compare()
+    results = load_data()
+    results_problems, results_types, results_master, results_nightly = clean_data(results)
+    table = Matrix{Any}(undef, length(results_master), 4)
+    fail = false
+    tolerance = 0.02
+    for (i, (master, nightly)) in enumerate(zip(results_master, results_nightly))
+        if results_types[i] == "time"
+            master = 1e9 .* master
+            nightly = 1e9 .* nightly
+            f, unit = best_unit(maximum(master))
+            m1 = round(mean(master) / f, digits = 2)
+            d1 = round(std(master) / f, digits = 2)
+            label_master = "$m1 ± $d1 $unit"
+            m2 = round(mean(nightly) / f, digits = 2)
+            d2 = round(std(nightly) / f, digits = 2)
+            label_nightly = "$m2 ± $d2 $unit"
+            indicator = if mean(master) < 1e9 * IGNORE_SMALL
+                0, "insignificant"
+            elseif (1 + MAX_DEVIATION) * m1 < m2
+                fail = true
+                2, "worse❌"
+            elseif m1 > (1 + MAX_DEVIATION) * m2
+                1, "better✅"
+            else
+                0, "don't care"
+            end
+        elseif results_types[i] == "allocs"
+            label_master = mean(master)
+            label_nightly = mean(nightly)
+            indicator = if label_master < (1 - tolerance) * label_nightly
+                fail = true
+                2, "worse❌"
+            elseif label_master > (1 + tolerance) * label_nightly
+                1, "better✅"
+            else
+                0, "don't care"
+            end
+        else
+            error("Beda!")
+        end
+        table[i, 1] = results_problems[i]
+        table[i, 2] = label_master
+        table[i, 3] = label_nightly
+        table[i, 4] = indicator[2]
+    end
+    fail, table
+end
+
+function post(fail, table)
+    comment_header = """
+    ## Running times benchmark
+
+    Note, that these numbers may fluctuate on the CI servers, so take them with a grain of salt.
+
+    """
+    io = IOBuffer()
+    println(io, comment_header)
+    if fail
+        println(io, "Potential regressions detected❌")
+    else
+        println(io, "No regressions detected✅")
+    end
+    table_header = ["Problem", "Master", "This commit", "Result"]
+    pretty_table(io, table, header = table_header, alignment = [:l, :r, :r, :r])
+    comment_str = String(take!(io))
+    println(comment_str)
+end
+
+function main()
+    runbench()
+    fail, table = compare()
+    post(fail, table)
+    versioninfo(verbose = true)
+    @testset "Benchmarks" begin
+        @test !fail
+        if fail
+            exit(1)
+        end
+    end
+end
+
+main()