diff --git a/src/ScalarEq/rhs.jl b/src/ScalarEq/rhs.jl
index 8f4f038abe838aee08e9f4da3d0be00f1815e35e..9cd9b0e5b2e4c512451566262cad9f68f9442d42 100644
--- a/src/ScalarEq/rhs.jl
+++ b/src/ScalarEq/rhs.jl
@@ -268,7 +268,7 @@ end
 
 
 function rhs!(env, mesh::Mesh1d, P::Project, hrsc::Nothing)
-  if mesh.element.kind === :modal_bspline2
+  if mesh.element.kind ∈ (:modal_bspline2, :modal_bspline1)
     modal_rhs!(env, mesh, P, hrsc)
   else
     nodal_rhs!(env, mesh, P, hrsc)
diff --git a/src/ScalarEq/setup.jl b/src/ScalarEq/setup.jl
index 3f906e69283abeb977c2f9e16b40d41a8a53ccfb..727b4242262ea5d758db63360ae1e5230e1031c5 100644
--- a/src/ScalarEq/setup.jl
+++ b/src/ScalarEq/setup.jl
@@ -14,8 +14,12 @@ function Project(env::Environment, mesh::Mesh1d, prms)
   analyze_error = prms["ScalarEq"]["analyze_error"]
   analyze_error_norm = prms["ScalarEq"]["analyze_error_norm"]
   V_bernstein = dg1d.Bernstein.vandermonde_matrix(mesh.element.z)
-  if mesh.element isa DGElement && mesh.element.kind === :modal_bspline2
-    V_bspline = dg1d.Bspline.vandermonde_matrix(mesh.element.z, mesh.element.bs2)
+  if mesh.element isa DGElement && mesh.element.kind ∈ (:modal_bspline2, :modal_bspline1)
+    if mesh.element.kind === :modal_bspline2
+      V_bspline = dg1d.Bspline.vandermonde_matrix(mesh.element.z, mesh.element.bs2)
+    elseif mesh.element.kind === :modal_bspline1
+      V_bspline = dg1d.Bspline.vandermonde_matrix(mesh.element.z, mesh.element.bs1)
+    end
     invV_bspline = pinv(V_bspline)
   else
     V_bspline, invV_bspline = nothing, nothing
@@ -132,7 +136,7 @@ function register_equation!(mesh::Mesh1d, P::Project)
       TODO(P.prms.hrsc)
     end
   end
-  if mesh.element isa DGElement && mesh.element.kind === :modal_bspline2
+  if mesh.element isa DGElement && mesh.element.kind ∈ (:modal_bspline2, :modal_bspline1)
     register_variables!(mesh.cache, static_variablenames=(:u_modal,:flx_u_modal,:src_u_modal))
   end
 end
diff --git a/src/bspline.jl b/src/bspline.jl
index 96e5dd1c6c20d1d737acd1a5f914a032734a97ed..ba72694af931c2054f470a29c9ce6a5bd254a81d 100644
--- a/src/bspline.jl
+++ b/src/bspline.jl
@@ -362,7 +362,7 @@ function integrate(f::Function, idxs::AbstractVector{<:Integer}, bs::Bspline1)
   @toggled_assert length(idxs) == bs.Npts
   sum = 0.0
   @inbounds for i in bs.Npts
-    fi, xl, xr = f[idxs[i]], bs.xs[i], bs.xs[i+2]
+    fi, xl, xr = f(idxs[i]), bs.xs[i], bs.xs[i+2]
     sum += fi * (xr-xl)/2
   end
   return sum
diff --git a/src/dgelement.jl b/src/dgelement.jl
index 3cd707001f81dbe0aed61b84da9cffac751cd753..e78129c0af170cd6da1a62171086a2f493defe3a 100644
--- a/src/dgelement.jl
+++ b/src/dgelement.jl
@@ -83,7 +83,7 @@ struct DGElement
       purge_zeros!(quadr_z)
       z = collect(range(-1.0,1.0,Npts))
       purge_zeros!(z)
-      D = [ Bernstein.derivative_bernstein_polynomial(z[i], j-1, N) for i in 1:Npts, j in 1:Npts ]
+      D = [ Bernstein.derivative_bernstein_polynomial(z[i], j, N) for i in 1:Npts, j in 1:Npts ]
       M = Bernstein.mass_matrix(Npts)
       S = Bernstein.stiffness_matrix(Npts)
       invM = inv(M)
@@ -110,6 +110,7 @@ struct DGElement
       # for which a 4th order accurate LG rule should suffice.
       # We use a Gauss rule to avoid evaluating the Bspline on subcell interfaces and bdrys.
       quadr_z, quadr_w = LG.rule(3)
+      # TODO Derivatives at bspline control points is ill-defined, because the splines are C^0 there
       D = [ Bspline.derivative_polynomial(z[i], j, bs1) for i in 1:Npts, j in 1:Npts ]
       M = Bspline.mass_matrix(bs1)
       S = Bspline.stiffness_matrix(bs1)
@@ -210,8 +211,10 @@ end
 
 
 @inline function integrate(u::AbstractVector, el::DGElement)
-  if el.kind ∈ (:modal_bspline2, :modal_bspline1)
-    return Bspline.integrate(u, el.bs)
+  if el.kind === :modal_bspline1
+    return Bspline.integrate(u, el.bs1)
+  elseif el.kind === :modal_bspline2
+    return Bspline.integrate(u, el.bs2)
   else
     if el.quadr === :lgl || el.quadr === :lgl_lumped
       return LGL.integrate(u, el.quadr_z, el.quadr_w)
@@ -223,8 +226,10 @@ end
 
 
 @inline function integrate(fn::Function, idxs::AbstractVector{<:Integer}, el::DGElement)
-  if el.kind ∈ (:modal_bspline2, :modal_bspline1)
-    return Bspline.integrate(fn, idxs, el.bs)
+  if el.kind === :modal_bspline1
+    return Bspline.integrate(fn, idxs, el.bs1)
+  elseif el.kind === :modal_bspline2
+    return Bspline.integrate(fn, idxs, el.bs2)
   else
     if el.quadr === :lgl || el.quadr === :lgl_lumped
       return LGL.integrate(fn, idxs, el.quadr_z, el.quadr_w)
diff --git a/test/IntegrationTests/refs/modal_bspline1_advection_sine/modal_bspline1_advection_sine.toml b/test/IntegrationTests/refs/modal_bspline1_advection_sine/modal_bspline1_advection_sine.toml
new file mode 100644
index 0000000000000000000000000000000000000000..8643fb4e18ebd51aa37adc8998c8e5ae3be45960
--- /dev/null
+++ b/test/IntegrationTests/refs/modal_bspline1_advection_sine/modal_bspline1_advection_sine.toml
@@ -0,0 +1,21 @@
+[Evolution]
+cfl = 0.5
+tend = 2.2
+
+[ScalarEq]
+analyze_error = true
+equation = "advection"
+id = "sine"
+
+[Output]
+substep_every_iteration = 200
+substep_variables = ["u", "u_err"]
+aligned_ts = "$(collect(range(0.1,2.2,step=0.1)))"
+variables1d = ["u", "u_err"]
+
+[Mesh]
+dg_kind = "modal_bspline1"
+range = [0.0, 1.0]
+k = 30
+basis = "lgl"
+n = 5
diff --git a/test/IntegrationTests/refs/modal_bspline1_advection_sine/output1d.h5 b/test/IntegrationTests/refs/modal_bspline1_advection_sine/output1d.h5
new file mode 100644
index 0000000000000000000000000000000000000000..94502a788fa15853819625290ce8e69faf1860c2
Binary files /dev/null and b/test/IntegrationTests/refs/modal_bspline1_advection_sine/output1d.h5 differ
diff --git a/test/IntegrationTests/refs/modal_bspline1_advection_sine/substep_output.h5 b/test/IntegrationTests/refs/modal_bspline1_advection_sine/substep_output.h5
new file mode 100644
index 0000000000000000000000000000000000000000..98c378fc817b9d2abb2b7ae8c42b1a0e39c1eac5
Binary files /dev/null and b/test/IntegrationTests/refs/modal_bspline1_advection_sine/substep_output.h5 differ
diff --git a/test/IntegrationTests/refs/testconfig.toml b/test/IntegrationTests/refs/testconfig.toml
index f9687bb2450cee77da23cdc94519931fcb9c79f3..fd03eac817b42f2d7a39325bf49540e2cf9e6b91 100644
--- a/test/IntegrationTests/refs/testconfig.toml
+++ b/test/IntegrationTests/refs/testconfig.toml
@@ -10,6 +10,9 @@ variables1d = [ "u", "u_err" ]
 [modal_bernstein_advection_sine]
 variables1d = [ "u" ]
 
+[modal_bspline1_advection_sine]
+variables1d = [ "u", "u_err" ]
+
 [modal_bspline2_advection_sine]
 variables1d = [ "u", "u_err" ]
 
diff --git a/test/IntegrationTests/src/IntegrationTests.jl b/test/IntegrationTests/src/IntegrationTests.jl
index 53bf4029c1b56ab3eb7950cca0fb9f813f7ba16e..84f8fbd6cd8290d97636e3673aa07a5c13a3b103 100644
--- a/test/IntegrationTests/src/IntegrationTests.jl
+++ b/test/IntegrationTests/src/IntegrationTests.jl
@@ -529,7 +529,9 @@ function impl_runtests(tests::Vector{String}=String[])
         (data, i, j) -> (j == 3 || j == 4) && !data[i,j],
         crayon"red"
     )
-    pretty_table(logs, hcat(1:length(tests), tests, isnothing.(runs_status), comparsion_success);
+    test_ran = isnothing.(runs_status)
+    test_passed = map(s -> ismissing(s) ? false : s, comparsion_success)
+    pretty_table(logs, hcat(1:length(tests), tests, test_ran, test_passed);
                  tf = tf_unicode_rounded,
                  alignment = [ :l, :l, :r, :r ],
                  highlighters = (hl_success,hl_fail),