diff --git a/dumux/experimental/timestepping/newmarkbeta.hh b/dumux/experimental/timestepping/newmarkbeta.hh
new file mode 100644
index 0000000000000000000000000000000000000000..507963d442eb54945326a6ec628c796f4a9c9f82
--- /dev/null
+++ b/dumux/experimental/timestepping/newmarkbeta.hh
@@ -0,0 +1,131 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+//
+// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
+// SPDX-License-Identifier: GPL-3.0-or-later
+//
+/*!
+ * \file
+ * \ingroup Experimental
+ * \brief Newmark-beta time integration scheme
+ */
+
+#ifndef DUMUX_TIMESTEPPING_NEWMARK_BETA_HH
+#define DUMUX_TIMESTEPPING_NEWMARK_BETA_HH
+
+#include <dumux/common/parameters.hh>
+
+namespace Dumux::Experimental {
+
+/*!
+ * \brief Newmark-beta time integration scheme
+ * \tparam Scalar The scalar type used for the solution
+ * \tparam Solution The solution type used for the variables
+ * \note Newmark (1959) "A method of computation for structural dynamics"
+ * Journal of the Engineering Mechanics Division,
+ * https://doi.org/10.1061/JMCEA3.0000098
+ * \note This is typically used for PDEs of the form M*a + C*v + K*x = f
+ * with M the mass matrix, C the damping matrix, K the stiffness matrix,
+ * x the displacement, v the velocity and a the acceleration.
+ * An example is dynamic mechanics.
+ */
+template<class Scalar, class Solution>
+class NewmarkBeta
+{
+ using Block = typename Solution::block_type;
+public:
+ //! Construct the Newmark-beta time integration scheme
+ //! with the given β and γ parameters
+ NewmarkBeta(const Scalar beta, const Scalar gamma)
+ : beta_(beta)
+ , gamma_(gamma)
+ {}
+
+ //! Construct the Newmark-beta time integration scheme
+ //! where parameters are read from the parameter file
+ NewmarkBeta()
+ : NewmarkBeta(
+ getParam<Scalar>("NewmarkBeta.Beta", 0.25),
+ getParam<Scalar>("NewmarkBeta.Gamma", 0.5)
+ ) {}
+
+ //! Initialize the time integration scheme with the current solution
+ //! while setting the velocity and acceleration to zero
+ void initialize(const Solution& x)
+ {
+ x_ = x;
+
+ // make sure the size is correct, but the values are zero
+ v_ = x; v_ = 0.0;
+ a_ = x; a_ = 0.0;
+ }
+
+ //! Initialize the time integration scheme with the current solution
+ void initialize(const Solution& x, const Solution& v, const Solution& a)
+ {
+ x_ = x;
+ v_ = v;
+ a_ = a;
+ }
+
+ //! Update with new position x
+ void update(const Scalar dt, const Solution& x)
+ {
+ for (std::size_t i = 0; i < x_.size(); ++i)
+ {
+ const auto xNew = x[i];
+ const auto aNew = acceleration(i, dt, xNew);
+ const auto vNew = velocity(i, dt, xNew, aNew);
+
+ x_[i] = xNew;
+ v_[i] = vNew;
+ a_[i] = aNew;
+ }
+ }
+
+ //! new a in terms of the old x, v, a, the new x and the time step size
+ Block acceleration(const std::size_t index, const Scalar dt, const Block& x) const
+ {
+ const auto& x0 = x_[index];
+ const auto& v0 = v_[index];
+ const auto& a0 = a_[index];
+
+ return (x - x0 - dt*v0)/(beta_*dt*dt) + (beta_ - 0.5)/beta_*a0;
+ }
+
+ //! new v in terms of the old v, a, the new x, a and the time step size
+ Block velocity(const std::size_t index, const Scalar dt, const Block& x, const Block& a) const
+ {
+ const auto& v0 = v_[index];
+ const auto& a0 = a_[index];
+
+ return v0 + dt*((1.0-gamma_)*a0 + gamma_*a);
+ }
+
+ //! new v in terms of the old v, a, the new x and the time step size
+ Block velocity(const std::size_t index, const Scalar dt, const Block& x) const
+ {
+ const auto a = acceleration(index, dt, x);
+ return velocity(index, dt, x, a);
+ }
+
+ //! current position
+ Scalar position(const std::size_t index) const
+ { return x_[index]; }
+
+ //! current velocity
+ Scalar velocity(const std::size_t index) const
+ { return v_[index]; }
+
+ //! current acceleration
+ Scalar acceleration(const std::size_t index) const
+ { return a_[index]; }
+
+private:
+ Scalar beta_, gamma_; //!< Newmark-beta parameters
+ Solution x_, v_, a_; //!< current position, velocity and acceleration
+};
+
+} // end namespace Dumux::Experimental
+
+#endif
diff --git a/test/experimental/timestepping/CMakeLists.txt b/test/experimental/timestepping/CMakeLists.txt
index ea918317f3cfd30b56d49fd6b1907223d04722b9..d4da32cf5e4f17cdfbca62f86f9c5ac40b798cd2 100644
--- a/test/experimental/timestepping/CMakeLists.txt
+++ b/test/experimental/timestepping/CMakeLists.txt
@@ -2,3 +2,4 @@
# SPDX-License-Identifier: GPL-3.0-or-later
dumux_add_test(SOURCES test_timestepmethods.cc LABELS unit timestepping experimental)
+dumux_add_test(SOURCES test_newmarkbeta.cc LABELS unit timestepping experimental)
diff --git a/test/experimental/timestepping/test_newmarkbeta.cc b/test/experimental/timestepping/test_newmarkbeta.cc
new file mode 100644
index 0000000000000000000000000000000000000000..ebbf49d71d418594c671d03bb9ac667cf41eb81f
--- /dev/null
+++ b/test/experimental/timestepping/test_newmarkbeta.cc
@@ -0,0 +1,86 @@
+//
+// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
+// SPDX-License-Identifier: GPL-3.0-or-later
+//
+#include <config.h>
+
+#include <iostream>
+#include <vector>
+
+#include <dune/istl/bvector.hh>
+#include <dune/common/exceptions.hh>
+
+#include <dumux/common/parameters.hh>
+#include <dumux/common/initialize.hh>
+#include <dumux/nonlinear/findscalarroot.hh>
+#include <dumux/io/gnuplotinterface.hh>
+
+#include <dumux/experimental/timestepping/newmarkbeta.hh>
+
+int main(int argc, char* argv[])
+{
+ using namespace Dumux;
+
+ // maybe initialize MPI and/or multithreading backend
+ Dumux::initialize(argc, argv);
+
+ // integrate d^2x/dt^2 = -x with Newmark-beta scheme
+ Experimental::NewmarkBeta<double, Dune::BlockVector<double>> newmark;
+ Dune::BlockVector<double> x(1); x = -1.0;
+ Dune::BlockVector<double> v(1); v = 0.0;
+ Dune::BlockVector<double> a(1); a = 1.0;
+ newmark.initialize(x, v, a);
+
+ const double dt = 0.1;
+ const double tEnd = 16*M_PI;
+ const int nSteps = tEnd/dt;
+
+ const bool showPlot = getParam<bool>("Plot", false);
+ Dumux::GnuplotInterface<double> plotter(showPlot);
+ std::vector<double> tValues, xValues, xExact, vValues, vExact, aValues, aExact;
+ tValues.reserve(nSteps);
+ xValues.reserve(nSteps); xExact.reserve(nSteps);
+ vValues.reserve(nSteps); vExact.reserve(nSteps);
+ aValues.reserve(nSteps); aExact.reserve(nSteps);
+
+ for (int i = 0; i < nSteps; ++i)
+ {
+ const double t = dt*i;
+
+ // mimic a nonlinear solver by finding the root of the residual (d^2x/dt^2 + x = 0)
+ const auto residual = [&](const auto x){ return newmark.acceleration(0, dt, x) + x; };
+ const auto xNew = findScalarRootNewton(x[0], residual);
+
+ x[0] = xNew;
+ // update with the new position
+ // after the update we have the current position, velocity and acceleration
+ newmark.update(dt, x);
+
+ tValues.push_back(t);
+ xValues.push_back(newmark.position(0));
+ xExact.push_back(-std::cos(t));
+ vValues.push_back(newmark.velocity(0));
+ vExact.push_back(std::sin(t));
+ aValues.push_back(newmark.acceleration(0));
+ aExact.push_back(std::cos(t));
+
+ if (i % 50 == 0)
+ {
+ plotter.resetPlot();
+ plotter.addDataSetToPlot(tValues, xValues, "x.dat", "w l");
+ plotter.addDataSetToPlot(tValues, xExact, "xe.dat", "w l");
+ plotter.addDataSetToPlot(tValues, vValues, "v.dat", "w l");
+ plotter.addDataSetToPlot(tValues, vExact, "ve.dat", "w l");
+ plotter.addDataSetToPlot(tValues, aValues, "a.dat", "w l");
+ plotter.addDataSetToPlot(tValues, aExact, "ae.dat", "w l");
+ plotter.plot("result");
+ }
+ }
+
+ // compare the final position with the exact solution
+ const double error = std::abs(x[0] + std::cos(tEnd));
+ if (error > 0.006)
+ DUNE_THROW(Dune::Exception, "Integration error too large: " << error);
+
+ return 0;
+}