diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_viscous_shock_newton_krylov.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_viscous_shock_newton_krylov.jl
new file mode 100644
index 00000000000..7ea78642724
--- /dev/null
+++ b/examples/p4est_2d_dgsem/elixir_navierstokes_viscous_shock_newton_krylov.jl
@@ -0,0 +1,183 @@
+using Trixi
+
+using OrdinaryDiffEqSDIRK
+using LinearSolve # For Jacobian-free Newton-Krylov (GMRES) solver
+using ADTypes # For automatic differentiation via finite differences
+
+# This is the classic 1D viscous shock wave problem with analytical solution
+# for a special value of the Prandtl number.
+# The original references are:
+#
+# - R. Becker (1922)
+#   Stoßwelle und Detonation.
+#   [DOI: 10.1007/BF01329605](https://doi.org/10.1007/BF01329605)
+#
+#   English translations:
+#   Impact waves and detonation. Part I.
+#   https://ntrs.nasa.gov/api/citations/19930090862/downloads/19930090862.pdf
+#   Impact waves and detonation. Part II.
+#   https://ntrs.nasa.gov/api/citations/19930090863/downloads/19930090863.pdf
+#
+# - M. Morduchow, P. A. Libby (1949)
+#   On a Complete Solution of the One-Dimensional Flow Equations
+#   of a Viscous, Head-Conducting, Compressible Gas
+#   [DOI: 10.2514/8.11882](https://doi.org/10.2514/8.11882)
+#
+#
+# The particular problem considered here is described in
+# - L. G. Margolin, J. M. Reisner, P. M. Jordan (2017)
+#   Entropy in self-similar shock profiles
+#   [DOI: 10.1016/j.ijnonlinmec.2017.07.003](https://doi.org/10.1016/j.ijnonlinmec.2017.07.003)
+
+### Fixed parameters ###
+
+# Special value for which nonlinear solver can be omitted
+# Corresponds essentially to fixing the Mach number
+alpha = 0.5
+# We want kappa = cp * mu = mu_bar to ensure constant enthalpy
+prandtl_number() = 3 / 4
+
+### Free choices: ###
+gamma() = 5 / 3
+
+# In Margolin et al., the Navier-Stokes equations are given for an
+# isotropic stress tensor τ, i.e., ∇ ⋅ τ = μ Δu
+mu_isotropic() = 0.15
+mu_bar() = mu_isotropic() / (gamma() - 1) # Re-scaled viscosity
+
+rho_0() = 1
+v() = 1 # Shock speed
+
+domain_length = 4.0
+
+### Derived quantities ###
+
+Ma() = 2 / sqrt(3 - gamma()) # Mach number for alpha = 0.5
+c_0() = v() / Ma() # Speed of sound ahead of the shock
+
+# From constant enthalpy condition
+p_0() = c_0()^2 * rho_0() / gamma()
+
+l() = mu_bar() / (rho_0() * v()) * 2 * gamma() / (gamma() + 1) # Appropriate length scale
+
+"""
+    initial_condition_viscous_shock(x, t, equations)
+
+Classic 1D viscous shock wave problem with analytical solution
+for a special value of the Prandtl number.
+The version implemented here is described in
+- L. G. Margolin, J. M. Reisner, P. M. Jordan (2017)
+  Entropy in self-similar shock profiles
+  [DOI: 10.1016/j.ijnonlinmec.2017.07.003](https://doi.org/10.1016/j.ijnonlinmec.2017.07.003)
+"""
+function initial_condition_viscous_shock(x, t, equations)
+    y = x[1] - v() * t # Translated coordinate
+
+    # Coordinate transformation. See eq. (33) in Margolin et al. (2017)
+    chi = 2 * exp(y / (2 * l()))
+
+    w = 1 + 1 / (2 * chi^2) * (1 - sqrt(1 + 2 * chi^2))
+
+    rho = rho_0() / w
+    u = v() * (1 - w)
+    p = p_0() * 1 / w * (1 + (gamma() - 1) / 2 * Ma()^2 * (1 - w^2))
+
+    return prim2cons(SVector(rho, u, 0, p), equations)
+end
+initial_condition = initial_condition_viscous_shock
+
+###############################################################################
+# semidiscretization of the ideal compressible Navier-Stokes equations
+
+equations = CompressibleEulerEquations2D(gamma())
+
+# Trixi implements the stress tensor in deviatoric form, thus we need to
+# convert the "isotropic viscosity" to the "deviatoric viscosity"
+mu_deviatoric() = mu_bar() * 3 / 4
+equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu_deviatoric(),
+                                                          Prandtl = prandtl_number(),
+                                                          gradient_variables = GradientVariablesEntropy())
+
+solver = DGSEM(polydeg = 3, surface_flux = flux_hlle)
+
+coordinates_min = (-domain_length / 2, -domain_length / 2)
+coordinates_max = (domain_length / 2, domain_length / 2)
+
+trees_per_dimension = (12, 3)
+mesh = P4estMesh(trees_per_dimension,
+                 polydeg = 3, initial_refinement_level = 0,
+                 coordinates_min = coordinates_min, coordinates_max = coordinates_max,
+                 periodicity = (false, true))
+
+### Inviscid boundary conditions ###
+
+# Prescribe pure influx based on initial conditions
+function boundary_condition_inflow(u_inner, normal_direction::AbstractVector, x, t,
+                                   surface_flux_function,
+                                   equations::CompressibleEulerEquations2D)
+    u_cons = initial_condition_viscous_shock(x, t, equations)
+    flux = Trixi.flux(u_cons, normal_direction, equations)
+
+    return flux
+end
+
+boundary_conditions = Dict(:x_neg => boundary_condition_inflow,
+                           :x_pos => boundary_condition_do_nothing)
+
+### Viscous boundary conditions ###
+# For the viscous BCs, we use the known analytical solution
+velocity_bc = NoSlip() do x, t, equations_parabolic
+    Trixi.velocity(initial_condition_viscous_shock(x,
+                                                   t,
+                                                   equations_parabolic),
+                   equations_parabolic)
+end
+
+heat_bc = Isothermal() do x, t, equations_parabolic
+    Trixi.temperature(initial_condition_viscous_shock(x,
+                                                      t,
+                                                      equations_parabolic),
+                      equations_parabolic)
+end
+
+boundary_condition_parabolic = BoundaryConditionNavierStokesWall(velocity_bc, heat_bc)
+
+boundary_conditions_parabolic = Dict(:x_neg => boundary_condition_parabolic,
+                                     :x_pos => boundary_condition_parabolic)
+
+semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
+                                             initial_condition, solver;
+                                             boundary_conditions = (boundary_conditions,
+                                                                    boundary_conditions_parabolic))
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+alive_callback = AliveCallback(alive_interval = 1)
+analysis_callback = AnalysisCallback(semi, interval = 100)
+
+callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback)
+
+###############################################################################
+# run the simulation
+
+# Tolerances for GMRES residual, see https://jso.dev/Krylov.jl/stable/solvers/unsymmetric/#Krylov.gmres
+atol_lin_solve = 1e-5
+rtol_lin_solve = 1e-5
+
+# Jacobian-free Newton-Krylov (GMRES) solver
+linsolve = KrylovJL_GMRES(atol = atol_lin_solve, rtol = rtol_lin_solve)
+
+# Use (diagonally) implicit Runge-Kutta, see
+# https://docs.sciml.ai/DiffEqDocs/stable/tutorials/advanced_ode_example/#Using-Jacobian-Free-Newton-Krylov
+ode_alg = Kvaerno4(autodiff = AutoFiniteDiff(), linsolve = linsolve)
+
+atol_ode_solve = 1e-4
+rtol_ode_solve = 1e-4
+sol = solve(ode, ode_alg;
+            abstol = atol_ode_solve, reltol = rtol_ode_solve,
+            ode_default_options()..., callback = callbacks);
diff --git a/examples/tree_1d_dgsem/elixir_diffusion_ldg.jl b/examples/tree_1d_dgsem/elixir_diffusion_ldg.jl
index 68f5c650c01..143a9cc9840 100644
--- a/examples/tree_1d_dgsem/elixir_diffusion_ldg.jl
+++ b/examples/tree_1d_dgsem/elixir_diffusion_ldg.jl
@@ -44,15 +44,14 @@ boundary_conditions_parabolic = boundary_condition_periodic
 solver_parabolic = ViscousFormulationLocalDG()
 semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
                                              initial_condition,
-                                             solver;
-                                             solver_parabolic,
+                                             solver; solver_parabolic,
                                              boundary_conditions = (boundary_conditions,
                                                                     boundary_conditions_parabolic))
 
 ###############################################################################
 # ODE solvers, callbacks etc.
 
-# Create ODE problem with time span from 0.0 to 1.0
+# Create ODE problem with time span from 0.0 to 0.1
 tspan = (0.0, 0.1)
 ode = semidiscretize(semi, tspan)
 
diff --git a/examples/tree_1d_dgsem/elixir_diffusion_ldg_newton_krylov.jl b/examples/tree_1d_dgsem/elixir_diffusion_ldg_newton_krylov.jl
new file mode 100644
index 00000000000..8d20593bc0e
--- /dev/null
+++ b/examples/tree_1d_dgsem/elixir_diffusion_ldg_newton_krylov.jl
@@ -0,0 +1,71 @@
+using Trixi
+
+using OrdinaryDiffEqSDIRK
+using LinearSolve # For Jacobian-free Newton-Krylov (GMRES) solver
+using ADTypes # For automatic differentiation via finite differences
+
+###############################################################################
+# semidiscretization of the linear (advection) diffusion equation
+
+advection_velocity = 0.0 # Note: This renders the equation mathematically purely parabolic
+equations = LinearScalarAdvectionEquation1D(advection_velocity)
+diffusivity() = 0.5
+equations_parabolic = LaplaceDiffusion1D(diffusivity(), equations)
+
+# surface flux does not matter for pure diffusion problem
+solver = DGSEM(polydeg = 3, surface_flux = flux_central)
+
+coordinates_min = -convert(Float64, pi)
+coordinates_max = convert(Float64, pi)
+
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 4,
+                n_cells_max = 30_000)
+
+function initial_condition_pure_diffusion_1d_convergence_test(x, t,
+                                                              equation)
+    nu = diffusivity()
+    c = 0
+    A = 1
+    omega = 1
+    scalar = c + A * sin(omega * sum(x)) * exp(-nu * omega^2 * t)
+    return SVector(scalar)
+end
+initial_condition = initial_condition_pure_diffusion_1d_convergence_test
+
+solver_parabolic = ViscousFormulationLocalDG()
+semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
+                                             initial_condition,
+                                             solver; solver_parabolic)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 2.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+analysis_callback = AnalysisCallback(semi, interval = 10)
+alive_callback = AliveCallback(alive_interval = 1)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback)
+
+###############################################################################
+# run the simulation
+
+# Tolerances for GMRES residual, see https://jso.dev/Krylov.jl/stable/solvers/unsymmetric/#Krylov.gmres
+atol_lin_solve = 1e-6
+rtol_lin_solve = 1e-5
+
+# Jacobian-free Newton-Krylov (GMRES) solver
+linsolve = KrylovJL_GMRES(atol = atol_lin_solve, rtol = rtol_lin_solve)
+
+# Use (diagonally) implicit Runge-Kutta, see
+# https://docs.sciml.ai/DiffEqDocs/stable/tutorials/advanced_ode_example/#Using-Jacobian-Free-Newton-Krylov
+ode_alg = KenCarp47(autodiff = AutoFiniteDiff(), linsolve = linsolve)
+
+atol_ode_solve = 1e-5
+rtol_ode_solve = 1e-4
+sol = solve(ode, ode_alg;
+            abstol = atol_ode_solve, reltol = rtol_ode_solve,
+            ode_default_options()..., callback = callbacks);
diff --git a/test/Project.toml b/test/Project.toml
index c996ec756a5..b08dc1a8dfb 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -12,6 +12,7 @@ ECOS = "e2685f51-7e38-5353-a97d-a921fd2c8199"
 ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
 NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
 OrdinaryDiffEqFeagin = "101fe9f7-ebb6-4678-b671-3a81e7194747"
@@ -43,6 +44,7 @@ ECOS = "1.1.2"
 ExplicitImports = "1.0.1"
 ForwardDiff = "0.10.36, 1"
 LinearAlgebra = "1"
+LinearSolve = "2.36.1, 3"
 MPI = "0.20.6"
 NLsolve = "4.5.1"
 OrdinaryDiffEqFeagin = "1"
diff --git a/test/test_parabolic_1d.jl b/test/test_parabolic_1d.jl
index 8969b06eb13..41692015764 100644
--- a/test/test_parabolic_1d.jl
+++ b/test/test_parabolic_1d.jl
@@ -58,6 +58,20 @@ end
     end
 end
 
+@trixi_testset "TreeMesh1D: elixir_diffusion_ldg_newton_krylov.jl" begin
+    @test_trixi_include(joinpath(examples_dir(), "tree_1d_dgsem",
+                                 "elixir_diffusion_ldg_newton_krylov.jl"),
+                        l2=[4.2710445174631516e-6], linf=[2.28491835256861e-5])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    let
+        t = sol.t[end]
+        u_ode = sol.u[end]
+        du_ode = similar(u_ode)
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+    end
+end
+
 @trixi_testset "TreeMesh1D: elixir_advection_diffusion_restart.jl" begin
     @test_trixi_include(joinpath(examples_dir(), "tree_1d_dgsem",
                                  "elixir_advection_diffusion_restart.jl"),
diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl
index 41450ffea1f..72825a37c4a 100644
--- a/test/test_parabolic_2d.jl
+++ b/test/test_parabolic_2d.jl
@@ -866,6 +866,32 @@ end
     end
 end
 
+@trixi_testset "P4estMesh2D: elixir_navierstokes_viscous_shock_newton_krylov.jl" begin
+    @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem",
+                                 "elixir_navierstokes_viscous_shock_newton_krylov.jl"),
+                        tspan=(0.0, 0.1),
+                        l2=[
+                            3.468233560427797e-5,
+                            2.64864594855224e-5,
+                            7.879490760481979e-10,
+                            2.8748482665365446e-5
+                        ],
+                        linf=[
+                            0.00018754529350140103,
+                            0.00014045634087878067,
+                            9.043610782328732e-9,
+                            0.00014499382160382268
+                        ])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    let
+        t = sol.t[end]
+        u_ode = sol.u[end]
+        du_ode = similar(u_ode)
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+    end
+end
+
 @trixi_testset "elixir_navierstokes_SD7003airfoil.jl" begin
     @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem",
                                  "elixir_navierstokes_SD7003airfoil.jl"),