diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 1847d192..d1425bb2 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -197,13 +197,16 @@ jobs:
           if [ "${{ matrix.os }}" = "macos-latest" ]; then
             HDF5_LIBRARY=$(brew --prefix hdf5)/lib/libhdf5.dylib
             MPI_SWITCH=(--mpi-library "$(brew --prefix open-mpi)/lib/libmpi.dylib")
+            SO_EXT="dylib"
           else
             HDF5_LIBRARY=/usr/lib/x86_64-linux-gnu/hdf5/openmpi/libhdf5.so
             MPI_SWITCH=()
+            SO_EXT="so"
           fi
           mkdir libtrixi-julia
           cd libtrixi-julia
           cp ../install/share/libtrixi/LibTrixi.jl/examples/libelixir_* .
+          cp ../install/lib/libsource_terms."$SO_EXT" .
           mkdir julia-depot
           ../install/bin/libtrixi-init-julia .. \
               --hdf5-library "$HDF5_LIBRARY" \
@@ -284,6 +287,7 @@ jobs:
           # all controllers
           ../build/examples/trixi_controller_simple_c          . libelixir_tree1d_advection_basic.jl
           ../build/examples/trixi_controller_simple_f          . libelixir_tree1d_advection_basic.jl
+          ../build/examples/trixi_controller_simple_c          . libelixir_structured2d_source_terms_callback.jl
           ../build/examples/trixi_controller_data_load_c       . libelixir_t8code2d_euler_tracer_amr.jl
           ../build/examples/trixi_controller_data_load_f       . libelixir_t8code2d_euler_tracer_amr.jl
           ../build/examples/trixi_controller_data_store_c      . libelixir_t8code2d_euler_tracer_amr.jl
diff --git a/LibTrixi.jl/examples/libelixir_structured2d_source_terms_callback.jl b/LibTrixi.jl/examples/libelixir_structured2d_source_terms_callback.jl
new file mode 100644
index 00000000..e2bbf7c7
--- /dev/null
+++ b/LibTrixi.jl/examples/libelixir_structured2d_source_terms_callback.jl
@@ -0,0 +1,97 @@
+using LibTrixi
+using Trixi
+using OrdinaryDiffEqLowStorageRK
+using Libdl
+
+struct SourceTermsCallback
+    source_term_fptr::Ptr{Nothing}
+
+    function SourceTermsCallback(; so_path)
+        so_handle = dlopen(so_path)
+        @info "Opened library ", dlpath(so_handle)
+        source_term_fptr = dlsym(so_handle, "source_term_wave")
+        @info "Obtained function pointer ", source_term_fptr
+
+        new(source_term_fptr)
+    end
+end
+
+function (callback::SourceTermsCallback)(u, x, t, equations::CompressibleEulerEquations2D)
+    @unpack source_term_fptr = callback
+
+    n_vars = nvariables(equations)
+    buffer = Ref(ntuple(i -> 0.0, n_vars))
+    @ccall $source_term_fptr(u::Ptr{Cdouble}, x::Ptr{Cdouble}, t::Cdouble,
+                             equations.gamma::Cdouble, buffer::Ptr{Cdouble})::Cvoid
+    return SVector(buffer[])
+end
+
+# The function to create the simulation state needs to be named `init_simstate`
+function init_simstate()
+
+    ###############################################################################
+    # semidiscretization of the compressible Euler equations
+
+    # use external function in locally compiled shared object
+    source_terms_callback = SourceTermsCallback(so_path = joinpath(@__DIR__,
+                                                                   "libsource_terms.so"))
+
+    equations = CompressibleEulerEquations2D(1.4)
+
+    initial_condition = initial_condition_convergence_test
+
+    solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
+
+    coordinates_min = (0.0, 0.0)
+    coordinates_max = (2.0, 2.0)
+
+    cells_per_dimension = (16, 16)
+
+    mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max,
+                          periodicity = true)
+
+    semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver;
+                                        source_terms = source_terms_callback,
+                                        boundary_conditions = boundary_condition_periodic)
+
+    ###############################################################################
+    # ODE solvers, callbacks etc.
+
+    tspan = (0.0, 2.0)
+    ode = semidiscretize(semi, tspan)
+
+    summary_callback = SummaryCallback()
+
+    analysis_interval = 100
+    analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+    alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+    save_solution = SaveSolutionCallback(interval = 100,
+                                        save_initial_solution = true,
+                                        save_final_solution = true,
+                                        solution_variables = cons2prim)
+
+    stepsize_callback = StepsizeCallback(cfl = 1.0)
+
+    callbacks = CallbackSet(summary_callback,
+                            analysis_callback, alive_callback,
+                            save_solution,
+                            stepsize_callback)
+
+    ###############################################################################
+    # create OrdinaryDiffEq's `integrator`
+
+    integrator = init(ode,
+                      CarpenterKennedy2N54(williamson_condition=false),
+                      dt=1.0, # will be overwritten by the stepsize_callback
+                      save_everystep=false,
+                      callback=callbacks);
+
+    ###############################################################################
+    # Create simulation state
+
+    simstate = SimulationState(semi, integrator)
+
+    return simstate
+end
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index be992e3c..835637a0 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -22,6 +22,13 @@ if (NOT LIBTRIXI_BUILD_FORTRAN_API)
     list( FILTER EXAMPLES EXCLUDE REGEX ".*\.f90" )
 endif()
 
+# exemplary external library providing function used to compute source terms
+add_library ( source_terms SHARED
+    source_terms.h
+    source_terms.c
+)
+
+install( TARGETS source_terms )
 foreach ( EXAMPLE ${EXAMPLES} )
 
     get_filename_component ( EXAMPLE_EXT ${EXAMPLE} EXT )
diff --git a/examples/source_terms.c b/examples/source_terms.c
new file mode 100644
index 00000000..3357b758
--- /dev/null
+++ b/examples/source_terms.c
@@ -0,0 +1,27 @@
+#include <stdio.h>
+#include <math.h>
+
+/* 
+ * Example for an external source term evaluation
+ * Diagonally running wave, taken from Trixi.jl's source_terms_convergence_test
+ */
+void source_term_wave(const double * u, const double * x, const double t,
+                      const double gamma, double * sourceterm) {
+
+    const double c = 2.0;
+    const double A = 0.1;
+    const double L = 2.0;
+    const double f = 1.0 / L;
+    const double omega = 2 * M_PI * f;
+
+    const double si = sin(omega * (x[0] + x[1] - t));
+    const double co = cos(omega * (x[0] + x[1] - t));
+    const double rho = c + A * si;
+    const double rho_x = omega * A * co;
+    const double tmp = (2 * rho - 1) * (gamma - 1);
+
+    sourceterm[0] = rho_x;
+    sourceterm[1] = rho_x * (1 + tmp);
+    sourceterm[2] = sourceterm[1];
+    sourceterm[3] = 2 * rho_x * (rho + tmp);
+}
diff --git a/examples/source_terms.h b/examples/source_terms.h
new file mode 100644
index 00000000..cdddc196
--- /dev/null
+++ b/examples/source_terms.h
@@ -0,0 +1,2 @@
+void source_term_wave(const double * u, const double * x, const double t,
+                      const double gamma, double * sourceterm);