diff --git a/Project.toml b/Project.toml
index 34581a22..fdd45303 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "QuantumCollocation"
 uuid = "0dc23a59-5ffb-49af-b6bd-932a8ae77adf"
 authors = ["Aaron Trowbridge <aaron.j.trowbridge@gmail.com> and contributors"]
-version = "0.3.2"
+version = "0.3.3"
 
 [deps]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
@@ -37,7 +37,7 @@ Interpolations = "0.15"
 Ipopt = "1.6"
 JLD2 = "0.5"
 MathOptInterface = "1.31"
-NamedTrajectories = "0.2"
+NamedTrajectories = "0.2.3"
 ProgressMeter = "1.10"
 Reexport = "1.2"
 Symbolics = "6.14"
diff --git a/docs/Project.toml b/docs/Project.toml
index c75e6d6a..3647548c 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -5,3 +5,6 @@ LiveServer = "16fef848-5104-11e9-1b77-fb7a48bbb589"
 NamedTrajectories = "538bc3a1-5ab9-4fc3-b776-35ca1e893e08"
 QuantumCollocation = "0dc23a59-5ffb-49af-b6bd-932a8ae77adf"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
+
+[compat]
+NamedTrajectories = "0.2.3"
\ No newline at end of file
diff --git a/docs/literate/man/ipopt_callbacks.jl b/docs/literate/man/ipopt_callbacks.jl
new file mode 100644
index 00000000..c69a1876
--- /dev/null
+++ b/docs/literate/man/ipopt_callbacks.jl
@@ -0,0 +1,88 @@
+# ```@meta
+# CollapsedDocStrings = true
+# ```
+# # IpOpt Callbacks
+
+# This page describes the callback functions that can be used with the IpOpt solver (in the future, may describe more general callback behavior).
+
+# ## Callbacks
+
+using QuantumCollocation
+using NamedTrajectories
+
+import ..QuantumStateSmoothPulseProblem
+import ..Callbacks
+
+# By default, IpOpt callbacks are called at each optimization step with the following signature:
+function full_argument_list_callback(
+    alg_mod::Cint,
+    iter_count::Cint,
+    obj_value::Float64,
+    inf_pr::Float64,
+    inf_du::Float64,
+    mu::Float64,
+    d_norm::Float64,
+    regularization_size::Float64,
+    alpha_du::Float64,
+    alpha_pr::Float64,
+    ls_trials::Cint,
+)
+    return true
+end
+
+# This gives the user access to some of the optimization state internals at each iteration.
+# A callback function with any subset of these arguments can be passed into the `solve!` function via the `callback` keyword argument see below.
+
+# The callback function can be used to stop the optimization early by returning `false`. The following callback when passed to `solve!` will stop the optimization after the first iteration:
+my_callback = (kwargs...) -> false
+
+# Single initial and target states
+# --------------------------------
+T = 50
+Δt = 0.2
+sys = QuantumSystem(0.1 * GATES[:Z], [GATES[:X], GATES[:Y]])
+ψ_init =  Vector{ComplexF64}([1.0, 0.0])
+ψ_target =  Vector{ComplexF64}([0.0, 1.0])
+
+prob = QuantumStateSmoothPulseProblem(
+    sys, ψ_init, ψ_target, T, Δt;
+    ipopt_options=IpoptOptions(print_level=1), 
+    piccolo_options=PiccoloOptions(verbose=false)
+)
+    
+
+# The callback function can be used to monitor the optimization progress, save intermediate results, or modify the optimization process.
+# For example, the following callback function saves the optimization trajectory at each iteration - this can be useful for debugging or plotting the optimization progress.
+# `trajectory_history_callback` from the `Callbacks` module
+callback, trajectory_history = QuantumCollocation.Callbacks.trajectory_history_callback(prob)
+solve!(prob, max_iter=20, callback=callback)
+
+# Save trajectory images into files which can be used to create a gif like the following: 
+for (iter, traj) in enumerate(trajectory_history)
+    str_index = lpad(iter, length(string(length(trajectory_history))), "0")
+    plot("./iteration-$str_index-trajectory.png", traj, [:ψ̃, :a], xlims=(-Δt, (T+5)*Δt), ylims=(ψ̃1 = (-2, 2), a = (-1.1, 1.1)))
+end
+
+# ![pulse optimization animation](../../assets/animation.gif)
+
+# Using a callback to get the best trajectory from all the optimization iterations
+sys2 = QuantumSystem(0.15 * GATES[:Z], [GATES[:X], GATES[:Y]])
+ψ_init2 =  Vector{ComplexF64}([0.0, 1.0])
+ψ_target2 =  Vector{ComplexF64}([1.0, 0.0])
+
+# Using other callbacks from the callback library 
+# --------------------------------
+# Callback used here is `best_rollout_fidelity_callback` which appends the best trajectories based on monotonically increasing fidelity of the rollout
+prob2 = QuantumStateSmoothPulseProblem(
+    sys2, ψ_init2, ψ_target2, T, Δt;
+    ipopt_options=IpoptOptions(print_level=1), 
+    piccolo_options=PiccoloOptions(verbose=false)
+)
+
+best_trajectory_callback, best_trajectory_list = best_rollout_fidelity_callback(prob2)
+solve!(prob2, max_iter=20, callback=best_trajectory_callback)
+# fidelity of the last iterate
+@show Losses.fidelity(prob2)
+
+# fidelity of the best iterate
+@show QuantumCollocation.fidelity(best_trajectory_list[end], prob2.system)
diff --git a/docs/make.jl b/docs/make.jl
index d56a1f17..fcf4d946 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -15,6 +15,7 @@ pages = [
     "Manual" => [
         "Problem Templates" => "generated/man/problem_templates.md",
         "Embedded Operators" => "generated/man/embedded_operators.md",
+        "Callbacks" => "generated/man/ipopt_callbacks.md",
     ],
     "Examples" => [
         "Two Qubit Gates" => "generated/examples/two_qubit_gates.md",
diff --git a/docs/src/assets/animation.gif b/docs/src/assets/animation.gif
new file mode 100644
index 00000000..2a79f44d
Binary files /dev/null and b/docs/src/assets/animation.gif differ
diff --git a/docs/src/assets/min_time_smooth_pulse_animation.gif b/docs/src/assets/min_time_smooth_pulse_animation.gif
new file mode 100644
index 00000000..97eb0654
Binary files /dev/null and b/docs/src/assets/min_time_smooth_pulse_animation.gif differ
diff --git a/docs/src/assets/pulse-optimization.gif b/docs/src/assets/pulse-optimization.gif
new file mode 100644
index 00000000..69e1adf7
Binary files /dev/null and b/docs/src/assets/pulse-optimization.gif differ
diff --git a/docs/src/assets/smooth_pulse_animation.gif b/docs/src/assets/smooth_pulse_animation.gif
new file mode 100644
index 00000000..d955946b
Binary files /dev/null and b/docs/src/assets/smooth_pulse_animation.gif differ
diff --git a/docs/src/generated/man/ipopt_callbacks.md b/docs/src/generated/man/ipopt_callbacks.md
new file mode 100644
index 00000000..04815ca3
--- /dev/null
+++ b/docs/src/generated/man/ipopt_callbacks.md
@@ -0,0 +1,126 @@
+```@meta
+EditURL = "../../../literate/man/ipopt_callbacks.jl"
+```
+
+```@meta
+CollapsedDocStrings = true
+```
+# IpOpt Callbacks
+
+This page describes the callback functions that can be used with the IpOpt solver (in the future, may describe more general callback behavior).
+
+## Callbacks
+
+````@example ipopt_callbacks
+using QuantumCollocation
+using NamedTrajectories
+
+import ..QuantumStateSmoothPulseProblem
+import ..Callbacks
+````
+
+By default, IpOpt callbacks are called at each optimization step with the following signature:
+
+````@example ipopt_callbacks
+function full_argument_list_callback(
+    alg_mod::Cint,
+    iter_count::Cint,
+    obj_value::Float64,
+    inf_pr::Float64,
+    inf_du::Float64,
+    mu::Float64,
+    d_norm::Float64,
+    regularization_size::Float64,
+    alpha_du::Float64,
+    alpha_pr::Float64,
+    ls_trials::Cint,
+)
+    return true
+end
+````
+
+This gives the user access to some of the optimization state internals at each iteration.
+A callback function with any subset of these arguments can be passed into the `solve!` function via the `callback` keyword argument see below.
+
+The callback function can be used to stop the optimization early by returning `false`. The following callback when passed to `solve!` will stop the optimization after the first iteration:
+
+````@example ipopt_callbacks
+my_callback = (kwargs...) -> false
+````
+
+Single initial and target states
+--------------------------------
+
+````@example ipopt_callbacks
+T = 50
+Δt = 0.2
+sys = QuantumSystem(0.1 * GATES[:Z], [GATES[:X], GATES[:Y]])
+ψ_init =  Vector{ComplexF64}([1.0, 0.0])
+ψ_target =  Vector{ComplexF64}([0.0, 1.0])
+
+prob = QuantumStateSmoothPulseProblem(
+    sys, ψ_init, ψ_target, T, Δt;
+    ipopt_options=IpoptOptions(print_level=1),
+    piccolo_options=PiccoloOptions(verbose=false)
+)
+````
+
+The callback function can be used to monitor the optimization progress, save intermediate results, or modify the optimization process.
+For example, the following callback function saves the optimization trajectory at each iteration - this can be useful for debugging or plotting the optimization progress.
+`trajectory_history_callback` from the `Callbacks` module
+
+````@example ipopt_callbacks
+callback, trajectory_history = QuantumCollocation.Callbacks.trajectory_history_callback(prob)
+solve!(prob, max_iter=20, callback=callback)
+````
+
+Save trajectory images into files which can be used to create a gif like the following:
+
+````@example ipopt_callbacks
+for (iter, traj) in enumerate(trajectory_history)
+    str_index = lpad(iter, length(string(length(trajectory_history))), "0")
+    plot("./iteration-$str_index-trajectory.png", traj, [:ψ̃, :a], xlims=(-Δt, (T+5)*Δt), ylims=(ψ̃1 = (-2, 2), a = (-1.1, 1.1)))
+end
+````
+
+![pulse optimization animation](../../assets/animation.gif)
+
+Using a callback to get the best trajectory from all the optimization iterations
+
+````@example ipopt_callbacks
+sys2 = QuantumSystem(0.15 * GATES[:Z], [GATES[:X], GATES[:Y]])
+ψ_init2 =  Vector{ComplexF64}([0.0, 1.0])
+ψ_target2 =  Vector{ComplexF64}([1.0, 0.0])
+````
+
+Using other callbacks from the callback library
+--------------------------------
+Callback used here is `best_rollout_fidelity_callback` which appends the best trajectories based on monotonically increasing fidelity of the rollout
+
+````@example ipopt_callbacks
+prob2 = QuantumStateSmoothPulseProblem(
+    sys2, ψ_init2, ψ_target2, T, Δt;
+    ipopt_options=IpoptOptions(print_level=1),
+    piccolo_options=PiccoloOptions(verbose=false)
+)
+
+best_trajectory_callback, best_trajectory_list = best_rollout_fidelity_callback(prob2)
+solve!(prob2, max_iter=20, callback=best_trajectory_callback)
+````
+
+fidelity of the last iterate
+
+````@example ipopt_callbacks
+@show Losses.fidelity(prob2)
+````
+
+fidelity of the best iterate
+
+````@example ipopt_callbacks
+@show QuantumCollocation.fidelity(best_trajectory_list[end], prob2.system)
+````
+
+---
+
+*This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*
+
diff --git a/src/QuantumCollocation.jl b/src/QuantumCollocation.jl
index 83b2e042..415f3d2c 100644
--- a/src/QuantumCollocation.jl
+++ b/src/QuantumCollocation.jl
@@ -72,5 +72,8 @@ include("problem_solvers.jl")
 include("plotting.jl")
 @reexport using .Plotting
 
+include("callbacks.jl")
+@reexport using .Callbacks
+
 
 end
diff --git a/src/callbacks.jl b/src/callbacks.jl
new file mode 100644
index 00000000..8e237a48
--- /dev/null
+++ b/src/callbacks.jl
@@ -0,0 +1,173 @@
+module Callbacks
+
+export best_rollout_fidelity_callback
+export best_unitary_rollout_fidelity_callback
+export trajectory_history_callback
+
+using NamedTrajectories
+using TestItemRunner
+
+using ..Losses
+using ..Problems: QuantumControlProblem, get_datavec
+using ..Rollouts
+
+
+function best_rollout_callback(prob::QuantumControlProblem, rollout::Function)
+    best_value = 0.0
+    best_trajectories = []
+
+    function callback(args...)
+        traj = NamedTrajectory(get_datavec(prob), prob.trajectory)
+        value = rollout(traj, prob.system)
+        if value > best_value
+            best_value = value
+            push!(best_trajectories, traj)
+        end
+        return true
+    end
+
+    return callback, best_trajectories
+end
+
+function best_rollout_fidelity_callback(prob::QuantumControlProblem)
+    return best_rollout_callback(prob, fidelity)
+end
+
+function best_unitary_rollout_fidelity_callback(prob::QuantumControlProblem)
+    return best_rollout_callback(prob, unitary_fidelity)
+end
+
+function trajectory_history_callback(prob::QuantumControlProblem)
+    trajectory_history = []
+    function callback(args...)
+        push!(trajectory_history, NamedTrajectory(get_datavec(prob), prob.trajectory))
+        return true
+    end
+
+    return callback, trajectory_history
+end
+
+# ========================================================================== # 
+
+@testitem "Callback returns false early stops" begin
+    using MathOptInterface
+    const MOI = MathOptInterface
+    include("../test/test_utils.jl")
+
+    prob = smooth_quantum_state_problem()
+
+    my_callback = (kwargs...) -> false
+
+    initial = fidelity(prob)
+    solve!(prob, max_iter=20, callback=my_callback)
+    final = fidelity(prob)
+
+    # callback forces problem to exit early as per Ipopt documentation
+    @test MOI.get(prob.optimizer, MOI.TerminationStatus()) == MOI.INTERRUPTED
+    @test initial ≈ final atol=1e-2
+end
+
+
+@testitem "Callback can get internal history" begin
+    using MathOptInterface
+    using NamedTrajectories
+    const MOI = MathOptInterface
+    include("../test/test_utils.jl")
+
+    prob = smooth_quantum_state_problem()
+
+    callback, trajectory_history = trajectory_history_callback(prob)
+
+    solve!(prob, max_iter=20, callback=callback)
+    @test length(trajectory_history) == 21
+end
+
+@testitem "Callback can get best state trajectory" begin
+    using MathOptInterface
+    using NamedTrajectories
+    const MOI = MathOptInterface
+    include("../test/test_utils.jl")
+
+    prob, system = smooth_quantum_state_problem(return_system=true)
+
+    callback, best_trajs = best_rollout_fidelity_callback(prob)
+    @test length(best_trajs) == 0
+
+    # measure fidelity
+    before = fidelity(prob)
+    solve!(prob, max_iter=20, callback=callback)
+
+    # length must increase if iterations are made
+    @test length(best_trajs) > 0
+    @test best_trajs[end] isa NamedTrajectory
+    
+    # fidelity ranking
+    after = fidelity(prob)
+    best = fidelity(best_trajs[end], system)
+    
+    @test before < after
+    @test before < best
+    @test after ≤ best
+end
+
+@testitem "Callback can get best unitary trajectory" begin
+    using MathOptInterface
+    using NamedTrajectories
+    const MOI = MathOptInterface
+    include("../test/test_utils.jl")
+
+    prob, system = smooth_unitary_problem(return_system=true)
+
+    callback, best_trajs = best_unitary_rollout_fidelity_callback(prob)
+    @test length(best_trajs) == 0
+
+    # measure fidelity
+    before = unitary_fidelity(prob)
+    solve!(prob, max_iter=20, callback=callback)
+
+    # length must increase if iterations are made
+    @test length(best_trajs) > 0
+    @test best_trajs[end] isa NamedTrajectory
+    
+    # fidelity ranking
+    after = unitary_fidelity(prob)
+    best = unitary_fidelity(best_trajs[end], system)
+    
+    @test before < after
+    @test before < best
+    @test after ≤ best
+end
+
+@testitem "Callback with full parameter test" begin
+    using MathOptInterface
+    using NamedTrajectories
+    const MOI = MathOptInterface
+    include("../test/test_utils.jl")
+
+    prob = smooth_quantum_state_problem()
+
+    obj_vals = []
+    function get_history_callback(
+        alg_mod::Cint,
+        iter_count::Cint,
+        obj_value::Float64,
+        inf_pr::Float64,
+        inf_du::Float64,
+        mu::Float64,
+        d_norm::Float64,
+        regularization_size::Float64,
+        alpha_du::Float64,
+        alpha_pr::Float64,
+        ls_trials::Cint,
+    )
+        push!(obj_vals, obj_value)
+        return iter_count < 3
+    end
+
+    solve!(prob, max_iter=20, callback=get_history_callback)
+
+    @test MOI.get(prob.optimizer, MOI.TerminationStatus()) == MOI.INTERRUPTED
+    @test length(obj_vals) == 4   # problem init, iter 1, iter 2, iter 3 (terminate)
+end
+
+end
\ No newline at end of file
diff --git a/src/problem_solvers.jl b/src/problem_solvers.jl
index d1e8e928..2169d976 100644
--- a/src/problem_solvers.jl
+++ b/src/problem_solvers.jl
@@ -9,8 +9,35 @@ using ..Options
 
 using NamedTrajectories
 using MathOptInterface
+using Ipopt
 const MOI = MathOptInterface
 
+
+"""
+   solve!(prob::QuantumControlProblem;
+        init_traj=nothing,
+        save_path=nothing,
+        max_iter=prob.ipopt_options.max_iter,
+        linear_solver=prob.ipopt_options.linear_solver,
+        print_level=prob.ipopt_options.print_level,
+        remove_slack_variables=false,
+        callback=nothing
+        # state_type=:unitary,
+        # print_fidelity=false,
+    )
+
+    Call optimization solver to solve the quantum control problem with parameters and callbacks.
+
+# Arguments
+- `prob::QuantumControlProblem`: The quantum control problem to solve.
+- `init_traj::NamedTrajectory`: Initial guess for the control trajectory. If not provided, a random guess will be generated.
+- `save_path::String`: Path to save the problem after optimization.
+- `max_iter::Int`: Maximum number of iterations for the optimization solver.
+- `linear_solver::String`: Linear solver to use for the optimization solver (e.g., "mumps", "paradiso", etc).
+- `print_level::Int`: Verbosity level for the solver.
+- `remove_slack_variables::Bool`: Remove slack variables from the trajectory after optimization.
+- `callback::Function`: Callback function to call during optimization steps.
+"""
 function solve!(
     prob::QuantumControlProblem;
     init_traj=nothing,
@@ -19,6 +46,7 @@ function solve!(
     linear_solver::String=prob.ipopt_options.linear_solver,
     print_level::Int=prob.ipopt_options.print_level,
     remove_slack_variables::Bool=false,
+    callback=nothing
     # state_type::Symbol=:unitary,
     # print_fidelity::Bool=false,
 )
@@ -35,6 +63,10 @@ function solve!(
     else
         set_trajectory!(prob)
     end
+    
+    if !isnothing(callback)
+        MOI.set(prob.optimizer, Ipopt.CallbackFunction(), callback)
+    end
 
     # if print_fidelity
     #     if state_type == :ket
diff --git a/test/test_utils.jl b/test/test_utils.jl
index ce5d309d..f4c3b5e9 100644
--- a/test/test_utils.jl
+++ b/test/test_utils.jl
@@ -116,3 +116,38 @@ function named_trajectory_type_1(; free_time=false)
         goal=goal
     )
 end
+
+function smooth_quantum_state_problem(; return_system::Bool=false)
+    T = 50
+    Δt = 0.2
+    sys = QuantumSystem(0.1 * PAULIS[:Z], [PAULIS[:X], PAULIS[:Y]])
+    ψ_init = ComplexF64[1.0, 0.0]
+    ψ_target = ComplexF64[0.0, 1.0]
+    prob = QuantumStateSmoothPulseProblem(
+        sys, ψ_init, ψ_target, T, Δt;
+        ipopt_options=IpoptOptions(print_level=1), 
+        piccolo_options=PiccoloOptions(verbose=false)
+    )
+    if return_system
+        return prob, sys
+    else
+        return prob
+    end
+end
+
+function smooth_unitary_problem(; return_system::Bool=false)
+    T = 50
+    Δt = 0.2
+    sys = QuantumSystem(0.1 * PAULIS[:Z], [PAULIS[:X], PAULIS[:Y]])
+    U_goal = GATES[:H]
+    prob = UnitarySmoothPulseProblem(
+        sys, U_goal, T, Δt;
+        ipopt_options=IpoptOptions(print_level=1), 
+        piccolo_options=PiccoloOptions(verbose=false)
+    )
+    if return_system
+        return prob, sys
+    else
+        return prob
+    end
+end