diff --git a/.gitignore b/.gitignore
index 7d095569..9bb70b68 100644
--- a/.gitignore
+++ b/.gitignore
@@ -8,4 +8,7 @@ deps/petsc-*.tar.gz
 docs/build/
 docs/site/
 Manifest.toml
-wrapping/.CondaPkg
\ No newline at end of file
+wrapping/.CondaPkg
+**/.vscode
+**/run/
+
diff --git a/examples/ex51.jl b/examples/ex51.jl
new file mode 100644
index 00000000..1d605349
--- /dev/null
+++ b/examples/ex51.jl
@@ -0,0 +1,239 @@
+using PETSc, MPI, Printf
+
+# Small ODE to test TS accuracy.
+#
+# The ODE
+#                 u1_t = cos(t),
+#                 u2_t = sin(u2)
+# with analytical solution
+#                 u1(t) = sin(t),
+#                 u2(t) = 2 * atan(exp(t) * tan(0.5))
+# is used to test the accuracy of TS schemes.
+
+# PETSc callbacks take a raw `void*` context pointer. We keep the active
+# `petsclib` inside a mutable Julia object so we can pass a stable reference
+# through `pointer_from_objref()` and recover it inside the RHS routine
+# to access the library-specific types, functions, and constants.
+mutable struct Ex51Context{PetscLib <: PETSc.LibPETSc.PetscLibType}
+    petsclib::PetscLib
+end
+
+"""
+    ex51_rhs!
+
+Right-hand side routine passed to `TSSetRHSFunction`, corresponding to
+`RHSFunction` in PETSc's `ex51.c`.
+"""
+function ex51_rhs!(
+    ::PETSc.LibPETSc.CTS, # optional argument `TS ts` not needed here
+    t::Float64,
+    u_ptr::PETSc.LibPETSc.CVec,
+    f_ptr::PETSc.LibPETSc.CVec,
+    ctx_ptr::Ptr{Cvoid},
+)::PETSc.LibPETSc.PetscErrorCode
+    # This function is called through PETSc's C callback interface, so its argument
+    # types must match the low-level PETSc ABI. In particular, the state and RHS
+    # vectors arrive as raw `CVec` pointers. We immediately wrap them using
+    # PETSc.jl's higher-level `VecPtr(..., own = false)` helper so we can use
+    # Julia-friendly helpers such as `withlocalarray!` without taking ownership
+    # away from PETSc.
+    ctx = unsafe_pointer_to_objref(ctx_ptr)::Ex51Context
+    petsclib = ctx.petsclib
+    # `own = false` since memory is managed by PETSc internally
+    u = PETSc.VecPtr(petsclib, u_ptr, false)
+    f = PETSc.VecPtr(petsclib, f_ptr, false)
+
+    PETSc.withlocalarray!(
+        (u, f);
+        read = (true, false),
+        write = (false, true),
+    ) do u_array, f_array
+        f_array[1] = cos(t)
+        f_array[2] = sin(u_array[2])
+    end
+
+    return PETSc.LibPETSc.PetscErrorCode(0)
+end
+
+const EX51_RHS_FUNCTION_PTR = @cfunction(
+    ex51_rhs!,
+    PETSc.LibPETSc.PetscErrorCode,
+    (
+        PETSc.LibPETSc.CTS,
+        Float64,
+        PETSc.LibPETSc.CVec,
+        PETSc.LibPETSc.CVec,
+        Ptr{Cvoid},
+    ),
+)
+
+function set_initial_condition!(u::PETSc.LibPETSc.PetscVec)
+    PETSc.withlocalarray!(u; read = false, write = true) do u_array
+        u_array[1] = 0.0
+        u_array[2] = 1.0
+    end
+    return nothing
+end
+
+function exact_solution!(u::PETSc.LibPETSc.PetscVec, t::Real)
+    PETSc.withlocalarray!(u; read = false, write = true) do u_array
+        u_array[1] = sin(t)
+        u_array[2] = 2 * atan(exp(t) * tan(0.5))
+    end
+    return nothing
+end
+
+"""
+    solve_ex51(; kwargs...)
+
+Solve the small ODE from PETSc TS tutorial `ex51.c` using PETSc.jl's low-level
+TS wrappers.
+
+Keyword arguments:
+
+- `petsclib`: PETSc library instance. Defaults to `Float64`.
+- `final_time`: final integration time. Default `1.0`.
+- `dt`: initial time step. Default `0.25`.
+- `save_trajectory`: whether to call `TSSetSaveTrajectory`. Default `true`.
+- `finalize_petsc`: whether to finalize PETSc at the end when this function had
+  to initialize it. Default `false` so repeated calls in one Julia session work
+  reliably.
+- `options`: additional PETSc command-line options for this solve. By default
+  the script uses `ARGS`, so invocations like
+  `julia --project=. examples/ex51.jl -ts_type rk -ts_rk_type 5dp` work, and
+  programmatic calls can use values such as
+  `options = ["-ts_type", "rk", "-ts_rk_type", "3bs"]`.
+- `verbose`: print a short run summary.
+
+Returns a named tuple with the final time, error norm, and numerical solution.
+"""
+function solve_ex51(;
+    petsclib = PETSc.getlib(PetscScalar = Float64),
+    final_time::Real = 1.0,
+    dt::Real = 0.25,
+    save_trajectory::Bool = true,
+    finalize_petsc::Bool = false,
+    options = nothing,
+    verbose::Bool = true,
+)
+    options === nothing && (options = copy(String.(ARGS)))
+    parsed_options = PETSc.parse_options(options)
+
+    comm = MPI.COMM_WORLD
+    PetscScalar = petsclib.PetscScalar
+
+    # `did_initialize`: whether we initialized the library in this call
+    did_initialize = !PETSc.initialized(petsclib)
+    if did_initialize
+        PETSc.initialize(petsclib)
+    end
+
+    MPI.Comm_size(comm) == 1 || error("This example only supports sequential runs.")
+
+    # Keep these variables concretely typed across the whole function, even
+    # though the actual PETSc objects are created later inside the `try` block.
+    # The null-pointer placeholders are overwritten by `TSCreate`/`VecSeq`,
+    # and the `finally` block checks `ptr != C_NULL` before destroying them.
+    ts = PETSc.LibPETSc.TS(petsclib)
+    u = PETSc.LibPETSc.PetscVec(petsclib)
+    u_exact = PETSc.LibPETSc.PetscVec(petsclib)
+    current_time = petsclib.PetscReal(NaN)
+    error_norm = petsclib.PetscReal(NaN)
+    solution = PetscScalar[]
+    ctx = Ex51Context(petsclib)
+    petsc_options = PETSc.Options(petsclib; parsed_options...)
+    pushed_options = false
+
+    try
+        # Create timestepping solver context.
+        ts = PETSc.LibPETSc.TSCreate(petsclib, comm)
+        PETSc.LibPETSc.TSSetType(petsclib, ts, "rosw")
+        PETSc.LibPETSc.TSSetProblemType(petsclib, ts, PETSc.LibPETSc.TS_NONLINEAR)
+
+        # Set initial conditions.
+        u = PETSc.VecSeq(petsclib, 2)
+        set_initial_condition!(u)
+        PETSc.LibPETSc.TSSetSolution(petsclib, ts, u)
+
+        callback_ctx_ptr = pointer_from_objref(ctx)
+        GC.@preserve ctx begin
+            PETSc.LibPETSc.TSSetRHSFunction(
+                petsclib,
+                ts,
+                nothing,
+                EX51_RHS_FUNCTION_PTR,
+                callback_ctx_ptr,
+            )
+
+            # Configure solver options. As in the C example, adaptivity is
+            # forced to take constant time steps.
+            save_trajectory && PETSc.LibPETSc.TSSetSaveTrajectory(petsclib, ts)
+            PETSc.LibPETSc.TSSetMaxTime(petsclib, ts, petsclib.PetscReal(final_time))
+            PETSc.LibPETSc.TSSetExactFinalTime(
+                petsclib,
+                ts,
+                PETSc.LibPETSc.TS_EXACTFINALTIME_STEPOVER,
+            )
+            PETSc.LibPETSc.TSSetTimeStep(petsclib, ts, petsclib.PetscReal(dt))
+
+            adapt = PETSc.LibPETSc.TSGetAdapt(petsclib, ts)
+            PETSc.LibPETSc.TSAdaptSetType(petsclib, adapt, "none")
+
+            push!(petsc_options)
+            pushed_options = true
+            PETSc.LibPETSc.TSSetFromOptions(petsclib, ts)
+            pop!(petsc_options)
+            pushed_options = false
+            PETSc.LibPETSc.TSSolve(petsclib, ts, u)
+        end
+
+        # Compute the error against the analytical solution at the achieved
+        # final time, matching the PETSc example's `VecAYPX` + `VecNorm` path.
+        current_time = PETSc.LibPETSc.TSGetTime(petsclib, ts)
+        u_exact = similar(u)
+        exact_solution!(u_exact, current_time)
+
+        PETSc.LibPETSc.VecAYPX(petsclib, u_exact, PetscScalar(-1), u)
+        error_norm = PETSc.LibPETSc.VecNorm(petsclib, u_exact, PETSc.NORM_2)
+        solution = copy(u[:])
+
+        if verbose
+            if abs(current_time - final_time) > 100 * eps(petsclib.PetscReal)
+                @printf(
+                    "Note: prescribed final time %.16g differs from actual final time %.16g\n",
+                    final_time,
+                    current_time,
+                )
+            end
+
+            @printf("TS type: %s\n", PETSc.LibPETSc.TSGetType(petsclib, ts))
+            @printf("Final time: %.16g\n", current_time)
+            @printf("Error at final time: %.2E\n", error_norm)
+        end
+
+        return (final_time = current_time, error = error_norm, solution = solution)
+    finally
+        if pushed_options
+            pop!(petsc_options)
+        end
+        if petsc_options.ptr != C_NULL
+            PETSc.destroy(petsc_options)
+        end
+        if u_exact.ptr != C_NULL
+            PETSc.destroy(u_exact)
+        end
+        if u.ptr != C_NULL
+            PETSc.destroy(u)
+        end
+        if ts.ptr != C_NULL
+            PETSc.LibPETSc.TSDestroy(petsclib, ts)
+        end
+        if did_initialize && finalize_petsc
+            PETSc.finalize(petsclib)
+        end
+    end
+end
+
+if !isinteractive() && abspath(PROGRAM_FILE) == @__FILE__
+    solve_ex51(finalize_petsc = true)
+end
diff --git a/src/PETSc.jl b/src/PETSc.jl
index 7a13b0ce..53a8d1f9 100644
--- a/src/PETSc.jl
+++ b/src/PETSc.jl
@@ -30,6 +30,7 @@ include("init.jl")
 include("vec.jl")       
 include("mat.jl")          
 include("options.jl")
+include("ts.jl")
 include("ksp.jl")
 include("snes.jl")          
 include("dm.jl")          
diff --git a/src/init.jl b/src/init.jl
index dd775a25..0e80d379 100644
--- a/src/init.jl
+++ b/src/init.jl
@@ -303,4 +303,4 @@ function check_petsc_wrappers_version(petsclib=nothing)
     end
 
     return (wrappers_version = wrappers_version, installed_version = installed_version, match = match)
-end
\ No newline at end of file
+end
diff --git a/src/ts.jl b/src/ts.jl
new file mode 100644
index 00000000..16cd0019
--- /dev/null
+++ b/src/ts.jl
@@ -0,0 +1,106 @@
+"""
+    TSSetRHSFunction(petsclib, ts, r, fptr::Ptr{Cvoid}, ctx = C_NULL)
+
+Convenience overload for low-level TS RHS callbacks created with `@cfunction`.
+
+The generated bindings currently accept the PETSc function-wrapper type directly,
+while Julia's `@cfunction` returns a raw pointer. This overload bridges that
+gap so callback-based TS examples can use the low-level interface naturally.
+"""
+function LibPETSc.TSSetRHSFunction(
+    petsclib::LibPETSc.PetscLibType,
+    ts::LibPETSc.TS,
+    r::AbstractPetscVec,
+    fptr::Ptr{Cvoid},
+    ctx::Ptr{Cvoid} = C_NULL,
+) end
+
+LibPETSc.@for_petsc function LibPETSc.TSSetRHSFunction(
+    petsclib::$UnionPetscLib,
+    ts::LibPETSc.TS,
+    r::AbstractPetscVec{$PetscLib},
+    fptr::Ptr{Cvoid},
+    ctx::Ptr{Cvoid} = C_NULL,
+)
+    typed_fptr = Base.unsafe_convert(Ptr{LibPETSc.TSRHSFunctionFn}, fptr)
+    LibPETSc.@chk ccall(
+        (:TSSetRHSFunction, $petsc_library),
+        LibPETSc.PetscErrorCode,
+        (LibPETSc.CTS, LibPETSc.CVec, Ptr{LibPETSc.TSRHSFunctionFn}, Ptr{Cvoid}),
+        ts,
+        r,
+        typed_fptr,
+        ctx,
+    )
+    return nothing
+end
+
+function LibPETSc.TSSetRHSFunction(
+    petsclib::LibPETSc.PetscLibType,
+    ts::LibPETSc.TS,
+    ::Nothing,
+    fptr::Ptr{Cvoid},
+    ctx::Ptr{Cvoid} = C_NULL,
+) end
+
+LibPETSc.@for_petsc function LibPETSc.TSSetRHSFunction(
+    petsclib::$UnionPetscLib,
+    ts::LibPETSc.TS,
+    ::Nothing,
+    fptr::Ptr{Cvoid},
+    ctx::Ptr{Cvoid} = C_NULL,
+)
+    typed_fptr = Base.unsafe_convert(Ptr{LibPETSc.TSRHSFunctionFn}, fptr)
+    LibPETSc.@chk ccall(
+        (:TSSetRHSFunction, $petsc_library),
+        LibPETSc.PetscErrorCode,
+        (LibPETSc.CTS, LibPETSc.CVec, Ptr{LibPETSc.TSRHSFunctionFn}, Ptr{Cvoid}),
+        ts,
+        C_NULL,
+        typed_fptr,
+        ctx,
+    )
+    return nothing
+end
+
+"""
+    adapt = TSGetAdapt(petsclib, ts)
+
+Return the adaptive time-step controller attached to `ts`.
+"""
+function LibPETSc.TSGetAdapt(
+    petsclib::LibPETSc.PetscLibType,
+    ts::LibPETSc.TS,
+) end
+
+LibPETSc.@for_petsc function LibPETSc.TSGetAdapt(
+    petsclib::$UnionPetscLib,
+    ts::LibPETSc.TS,
+)
+    adapt_ref = Ref{LibPETSc.TSAdapt}()
+    LibPETSc.@chk ccall(
+        (:TSGetAdapt, $petsc_library),
+        LibPETSc.PetscErrorCode,
+        (LibPETSc.CTS, Ptr{LibPETSc.TSAdapt}),
+        ts,
+        adapt_ref,
+    )
+    return adapt_ref[]
+end
+
+"""
+    TSAdaptSetType(petsclib, adapt, type::String)
+
+Convenience wrapper for setting the TS adaptivity controller using a Julia
+string such as `"none"` or `"basic"`.
+"""
+function LibPETSc.TSAdaptSetType(
+    petsclib::LibPETSc.PetscLibType,
+    adapt::LibPETSc.TSAdapt,
+    type::String,
+)
+    c_str = Vector{UInt8}(type * "\0")
+    ptr = Base.unsafe_convert(Ptr{Cchar}, pointer(c_str))
+    LibPETSc.TSAdaptSetType(petsclib, adapt, ptr)
+    return nothing
+end
diff --git a/test/examples.jl b/test/examples.jl
index 87368712..666eba5f 100644
--- a/test/examples.jl
+++ b/test/examples.jl
@@ -5,7 +5,7 @@ using MPI
 using .PETScTestUtils: find_sources
 
 @testset "examples" begin
-  examples_dir = joinpath(@__DIR__, "..", "examples")
+  examples_dir = joinpath(dirname(@__DIR__), "examples")
   examples = find_sources(examples_dir)
   filter!(file -> readline(file) != "# EXCLUDE FROM TESTING", examples)
 
diff --git a/test/mpi_examples.jl b/test/mpi_examples.jl
index 18a838ce..7cc0216c 100644
--- a/test/mpi_examples.jl
+++ b/test/mpi_examples.jl
@@ -4,7 +4,7 @@ using MPI
 using .PETScTestUtils: find_sources
 
 @testset "mpi examples" begin
-  examples_dir = joinpath(@__DIR__, "..", "examples")
+  examples_dir = joinpath(dirname(@__DIR__), "examples")
   examples = find_sources(examples_dir)
   filter!(file -> readline(file) == "# INCLUDE IN MPI TEST", examples)
 
diff --git a/test/runtests.jl b/test/runtests.jl
index a6f54027..2462a15e 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -30,11 +30,12 @@ include("dmnetwork.jl")     # new test for DMNetwork example
 include("dmshell.jl")       # new test for DMShell example
 include("dmproduct.jl")     # test for DMProduct example
 include("matshell.jl")      # autowrapped!
-include("test_dmstag.jl") 
-include("test_snes.jl")  
+include("test_dmstag.jl")
+include("test_snes.jl")
 include("old_test.jl")
 include("low_level_viewer.jl")  # Low-level viewer convenience functions
 include("low_level_ts.jl")      # Low-level TS functions
+include("ts_ex51.jl")           # Regression test for repeated ex51 solves
 include("low_level_is.jl")      # Low-level IS functions
 include("low_level_petscsection.jl")  # Low-level PetscSection functions
 include("low_level_tao.jl")     # Low-level Tao functions
@@ -67,5 +68,3 @@ if do_mpi
         end
     end
 end
-
-
diff --git a/test/ts_ex51.jl b/test/ts_ex51.jl
new file mode 100644
index 00000000..28b8ca36
--- /dev/null
+++ b/test/ts_ex51.jl
@@ -0,0 +1,61 @@
+using Test
+using PETSc
+
+include(joinpath(dirname(@__DIR__), "examples", "ex51.jl"))
+
+@testset "TS ex51 example" begin
+    petsclib = PETSc.getlib(PetscScalar = Float64)
+    PETSc.initialize(petsclib)
+    try
+        # `examples.jl` covers the default `save_trajectory = true` path in a
+        # fresh Julia subprocess. Here we disable trajectory saving because
+        # PETSc 3.22 does not reliably re-register the built-in `basic`
+        # trajectory type after earlier TS usage and reinitialization in the
+        # same process.
+        result_default_1 = solve_ex51(;
+            petsclib,
+            options = String[],
+            save_trajectory = false,
+            verbose = false,
+        )
+        result_3bs = solve_ex51(;
+            petsclib,
+            options = ["-ts_type", "rk", "-ts_rk_type", "3bs"],
+            save_trajectory = false,
+            verbose = false,
+        )
+
+        result_default_2 = solve_ex51(;
+            petsclib,
+            options = String[],
+            save_trajectory = false,
+            verbose = false,
+        )
+        result_5dp = solve_ex51(;
+            petsclib,
+            options = ["-ts_type", "rk", "-ts_rk_type", "5dp"],
+            save_trajectory = false,
+            verbose = false,
+        )
+
+        @test result_default_1.final_time ≈ 1.0 atol = 100 * eps(Float64)
+        @test length(result_default_1.solution) == 2
+        @test result_default_1.error < 5.0e-4
+
+        @test result_3bs.final_time ≈ 1.0 atol = 100 * eps(Float64)
+        @test length(result_3bs.solution) == 2
+        @test result_3bs.error < 5.0e-5
+
+        @test result_default_2.final_time ≈ 1.0 atol = 100 * eps(Float64)
+        @test length(result_default_2.solution) == 2
+        @test result_default_2.error ≈ result_default_1.error rtol = 1e-12
+
+        @test result_5dp.final_time ≈ 1.0 atol = 100 * eps(Float64)
+        @test length(result_5dp.solution) == 2
+        @test result_5dp.error < 5.0e-8
+    finally
+        if PETSc.initialized(petsclib) && !PETSc.finalized(petsclib)
+            PETSc.finalize(petsclib)
+        end
+    end
+end