diff --git a/docs/conf.py b/docs/conf.py
index 439fb1b58..68ca060c2 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -30,13 +30,13 @@
     'sphinx.ext.mathjax',
     'sphinx.ext.viewcode',
     'sphinx.ext.extlinks',
-    'jupyter_sphinx',
     'sphinx_autodoc_typehints',
     'reno.sphinxext',
     'sphinx.ext.intersphinx',
     'nbsphinx',
     'sphinxcontrib.bibtex',
     "qiskit_sphinx_theme",
+    "matplotlib.sphinxext.plot_directive",
 ]
 templates_path = ["_templates"]
 
@@ -61,6 +61,10 @@
 
 html_context = {"version_list": ["0.4"]}
 
+html_theme_options = {
+    "sidebar_qiskit_ecosystem_member": True,
+}
+
 # autodoc/autosummary options
 autosummary_generate = True
 autosummary_generate_overwrite = False
diff --git a/docs/tutorials/Lindblad_dynamics_simulation.rst b/docs/tutorials/Lindblad_dynamics_simulation.rst
index 793706166..64c170706 100644
--- a/docs/tutorials/Lindblad_dynamics_simulation.rst
+++ b/docs/tutorials/Lindblad_dynamics_simulation.rst
@@ -64,7 +64,9 @@ syntax of ``Pauli`` classes to indicate a qubit number, as below.
 Below, we first set the number of qubits :math:`N` to be simulated, and then prepare and store the
 single-qubit Pauli operators that will be used in the rest of this tutorial.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     import numpy as np
     from qiskit.quantum_info import Operator, Pauli
@@ -100,7 +102,9 @@ two-qubit terms. Since there are no time-dependent terms, and we do not plan to
 derivatives of parameters, we do not use the :class:`Signal` class in this tutorial. See the other
 tutorials for various generalizations of this approach supported with ``qiskit-dynamics``.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     from qiskit_dynamics import Solver, Signal
 
@@ -140,7 +144,9 @@ tutorials for various generalizations of this approach supported with ``qiskit-d
 We now define the initial state for the simulation, the time span to simulate for, and the
 intermediate times for which the solution is requested.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     from qiskit.quantum_info import DensityMatrix
 
@@ -175,7 +181,9 @@ invariant as well. Hence the mean Bloch vector should be equal to any qubit’s
 observing that this equality holds is a simple and useful verification of the numerical solution
 that will be added in the next section.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     n_times = len(sol.y)
     x_data = np.zeros((N, n_times))
@@ -218,11 +226,12 @@ tilt along :math:`+x`, while for :math:`J=3` it will significantly shorten (the
 a mixed state), becoming tilted along :math:`-y`. This complex dependence of the Bloch vector on the
 parameters can be systematically analyzed - we encourage you to try it!
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     from qiskit.visualization import plot_bloch_vector
     import matplotlib.pyplot as plt
-    %matplotlib inline
 
     fontsize = 16
 
@@ -235,8 +244,8 @@ parameters can be systematically analyzed - we encourage you to try it!
     ax.set_xlabel('$t$', fontsize = fontsize)
     ax.set_title('Mean Bloch vector vs. $t$', fontsize = fontsize)
 
-    display(plot_bloch_vector([x_mean[-1], y_mean[-1], z_mean[-1]],
-                      f'Mean Bloch vector at $t = {t_eval[-1]}$'))
+    plot_bloch_vector([x_mean[-1], y_mean[-1], z_mean[-1]],
+                      f'Mean Bloch vector at $t = {t_eval[-1]}$')
 
     if N > 1 and ((abs(x_mean[-1]) > 1e-5 and abs(x_data[0, -1] / x_mean[-1] - 1) > 1e-5 or
                   (abs(z_mean[-1]) > 1e-5 and abs(z_data[1, -1] / z_mean[-1] - 1) > 1e-5))):
diff --git a/docs/tutorials/Rabi_oscillations.rst b/docs/tutorials/Rabi_oscillations.rst
index 729dfd757..ac69b0d4d 100644
--- a/docs/tutorials/Rabi_oscillations.rst
+++ b/docs/tutorials/Rabi_oscillations.rst
@@ -35,7 +35,9 @@ then setup the :class:`.Solver` class instance that stores and manipulates the m
 using matrices and :class:`.Signal` instances. For the time-independent :math:`z` term we set the
 signal to a constant, while for the trasverse driving term we setup a harmonic signal.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     import numpy as np
     from qiskit.quantum_info import Operator
@@ -62,7 +64,9 @@ signal to a constant, while for the trasverse driving term we setup a harmonic s
 We now define the initial state for the simulation, the time span to simulate for, and the
 intermediate times for which the solution is requested, and solve the evolution.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     from qiskit.quantum_info.states import Statevector
     from qiskit.quantum_info import DensityMatrix
@@ -96,11 +100,12 @@ increases and decreases). This mechanism of Rabi oscillations is the basis for t
 gates used to manipulate quantum devices - in particular this is a realization of the :math:`X`
 gate.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     from qiskit.visualization import plot_bloch_vector
     import matplotlib.pyplot as plt
-    %matplotlib inline
 
     fontsize = 16
 
@@ -125,8 +130,8 @@ gate.
         ax.set_title('Bloch vector vs. $t$', fontsize = fontsize)
         plt.show()
 
-        display(plot_bloch_vector([x_data[-1], y_data[-1], z_data[-1]],
-                                  f'Bloch vector at $t = {t_eval[-1]}$'))
+        plot_bloch_vector([x_data[-1], y_data[-1], z_data[-1]],
+                                  f'Bloch vector at $t = {t_eval[-1]}$')
 
     plot_qubit_dynamics(sol, t_eval, X, Y, Z)
 
@@ -161,7 +166,9 @@ since in ``qiskit`` and in ``qiskit-dynamics`` the syntax of many functions is i
 state vectors and density matrices. The shrinking of the qubit’s state within the Bloch sphere due
 to the incoherent evolution can be clearly seen in the plots below.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     Gamma_1 = .8
     Gamma_2 = .2
diff --git a/docs/tutorials/dynamics_backend.rst b/docs/tutorials/dynamics_backend.rst
index a78058248..fe9ae1e59 100644
--- a/docs/tutorials/dynamics_backend.rst
+++ b/docs/tutorials/dynamics_backend.rst
@@ -28,14 +28,16 @@ when using a JAX solver method. Here we configure JAX to run on CPU in 64 bit mo
 :ref:`User Guide entry on using different array libraries with Qiskit Dynamics <how-to use different
 array libraries>` for more information.
 
-.. jupyter-execute::
-    :hide-code:
+.. plot::
+    :context: close-figs
 
     # a parallelism warning raised by JAX is being raised due to somethign outside of Dynamics
     import warnings
     warnings.filterwarnings("ignore", message="os.fork")
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     # Configure JAX
     import jax
@@ -68,7 +70,9 @@ where
 - :math:`N_0` and :math:`N_1` are the number operators for qubits :math:`0` and :math:`1`
   respectively.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     import numpy as np
     
@@ -118,7 +122,9 @@ outcomes of :meth:`.DynamicsBackend.run` are independent of the choice of rotati
 :class:`.Solver`, and as such we are free to choose the rotating frame that provides the best
 performance.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     from qiskit_dynamics import Solver
     
@@ -146,7 +152,9 @@ Furthermore, note that in the solver options we set the max step size to the pul
 ``dt`` via the ``"hmax"`` argument for the method ``"jax_odeint"``. This is important for preventing
 variable step solvers from accidentally stepping over pulses in schedules with long idle times.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     from qiskit_dynamics import DynamicsBackend
     
@@ -182,12 +190,13 @@ that the usual instructions work on the :class:`.DynamicsBackend`.
     impact on the returned results.
 
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
-    %%time
-    
+    import time
     from qiskit import pulse
-    
+
     sigma = 128
     num_samples = 256
     
@@ -207,20 +216,28 @@ that the usual instructions work on the :class:`.DynamicsBackend`.
                 pulse.acquire(duration=1, qubit_or_channel=0, register=pulse.MemorySlot(0))
             
         schedules.append(schedule)
-        
+    
+    start_time = time.time()
+
     job = backend.run(schedules, shots=100)
     result = job.result()
 
+    print(f"Run time: {time.time() - start_time}")
+
 Visualize one of the schedules.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     schedules[3].draw()
 
 Retrieve the counts for one of the experiments as would be done using the results object from a real
 backend.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     result.get_counts(3)
 
@@ -238,7 +255,9 @@ contained in the :class:`.DynamicsBackend`.
 Build a simple circuit. Here we build one consisting of a single Hadamard gate on qubit :math:`0`,
 followed by measurement.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     from qiskit import QuantumCircuit
     
@@ -252,7 +271,9 @@ Next, attach a calibration for the Hadamard gate on qubit :math:`0` to the circu
 we are only demonstrating the mechanics of adding a calibration; we have not attempted to calibrate
 the schedule to implement the Hadamard gate with high fidelity.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     with pulse.build() as h_q0:
         pulse.play(
@@ -264,10 +285,16 @@ the schedule to implement the Hadamard gate with high fidelity.
 
 Call run on the circuit, and get counts as usual.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
-    %time res = backend.run(circ).result()
+    start_time = time.time()
     
+    res = backend.run(circ).result()
+    
+    print(f"Run time: {time.time() - start_time}")
+
     res.get_counts(0)
 
 4.2 Simulating circuits via gate definitions in the backend :class:`~qiskit.transpiler.Target`
@@ -277,7 +304,9 @@ Alternatively to the above work flow, add the above schedule as the pulse-level
 Hadamard gate on qubit :math:`0` to `backend.target`, which impacts how jobs are transpiled for the
 backend. See the :class:`~qiskit.transpiler.Target` class documentation for further information.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     from qiskit.circuit.library import HGate
     from qiskit.transpiler import InstructionProperties
@@ -287,15 +316,23 @@ backend. See the :class:`~qiskit.transpiler.Target` class documentation for furt
 Rebuild the same circuit, however this time we do not need to add the calibration for the Hadamard
 gate to the circuit object.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     circ2 = QuantumCircuit(1, 1)
     circ2.h(0)
     circ2.measure([0], [0])
     
-    %time result = backend.run(circ2).result()
+    start_time = time.time()
 
-.. jupyter-execute::
+    result = backend.run(circ2).result()
+
+    print(f"Run time: {time.time() - start_time}")
+
+.. plot::
+    :context: close-figs
+    :include-source:
 
     result.get_counts(0)
 
@@ -324,7 +361,9 @@ To enable running of the single qubit experiments, we add the following to the `
   be utilizing it, this ensures that validation steps checking that the device is fully connected 
   will pass.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     from qiskit.circuit.library import XGate, SXGate, RZGate, CXGate
     from qiskit.circuit import Parameter
@@ -364,7 +403,9 @@ object. Here we use the
 :class:`~qiskit_experiments.calibration_management.basis_gate_library.FixedFrequencyTransmon`
 template library to initialize our calibrations.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     import pandas as pd
     from qiskit_experiments.calibration_management.calibrations import Calibrations
@@ -380,7 +421,9 @@ template library to initialize our calibrations.
 Next, run a rough amplitude calibration for ``X`` and ``SX`` gates for both qubits. First, build the
 experiments.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     from qiskit_experiments.library.calibration import RoughXSXAmplitudeCal
     
@@ -392,25 +435,36 @@ experiments.
 
 Run the Rabi experiments.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
+
+    start_time = time.time()
 
-    %%time
     rabi0_data = rabi0.run().block_for_results()
     rabi1_data = rabi1.run().block_for_results()
 
+    print(f"Run time: {time.time() - start_time}")
+
 Plot the results.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     rabi0_data.figure(0)
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     rabi1_data.figure(0)
 
 Observe the updated parameters for qubit 0.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     pd.DataFrame(**cals.parameters_table(qubit_list=[0, ()], parameters="amp"))
 
@@ -420,7 +474,9 @@ Observe the updated parameters for qubit 0.
 Run rough Drag parameter calibration for the ``X`` and ``SX`` gates. This follows the same procedure
 as above.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     from qiskit_experiments.library.calibration import RoughDragCal
     
@@ -432,24 +488,35 @@ as above.
     
     cal_drag0.circuits()[5].draw(output="mpl")
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
+
+    start_time = time.time()
 
-    %%time
     drag0_data = cal_drag0.run().block_for_results()
     drag1_data = cal_drag1.run().block_for_results()
 
-.. jupyter-execute::
+    print(f"Run time: {time.time() - start_time}")
+
+.. plot::
+    :context: close-figs
+    :include-source:
 
     drag0_data.figure(0)
 
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     drag1_data.figure(0)
 
 The updated calibrations object:
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     pd.DataFrame(**cals.parameters_table(qubit_list=[0, ()], parameters="amp"))
 
@@ -463,7 +530,9 @@ channel map, which is a dictionary mapping control-target qubit index pairs (giv
 the control channel index used to drive the corresponding cross-resonance interaction. This is
 required by the experiment to determine which channel to drive for each control-target pair.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
     
     # set the control channel map
     backend.set_options(control_channel_map={(0, 1): 0, (1, 0): 1})
@@ -471,7 +540,9 @@ required by the experiment to determine which channel to drive for each control-
 Build the characterization experiment object, and update gate definitions in ``target`` with the
 values for the single qubit gates calibrated above.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     from qiskit_experiments.library import CrossResonanceHamiltonian
 
@@ -483,17 +554,27 @@ values for the single qubit gates calibrated above.
     
     backend.target.update_from_instruction_schedule_map(cals.get_inst_map())
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     cr_ham_experiment.circuits()[10].draw("mpl")
 
 Run the simulation.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
+
+    start_time = time.time()
+
+    data_cr = cr_ham_experiment.run().block_for_results()
 
-    %time data_cr = cr_ham_experiment.run().block_for_results()
+    print(f"Run time: {time.time() - start_time}")
 
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     data_cr.figure(0)
diff --git a/docs/tutorials/optimizing_pulse_sequence.rst b/docs/tutorials/optimizing_pulse_sequence.rst
index 2c3f25191..c681ddc4d 100644
--- a/docs/tutorials/optimizing_pulse_sequence.rst
+++ b/docs/tutorials/optimizing_pulse_sequence.rst
@@ -24,7 +24,9 @@ We will optimize an :math:`X`-gate on a model of a qubit system using the follow
 
 First, set JAX to operate in 64-bit mode and to run on CPU.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     import jax
     jax.config.update("jax_enable_x64", True)
@@ -51,7 +53,9 @@ In the above:
 
 We will setup the problem to be in the rotating frame of the drift term.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     import numpy as np
     from qiskit.quantum_info import Operator
@@ -100,7 +104,9 @@ We remark that there are many other parameterizations that may achieve the same
 more efficient strategies for achieving a value of :math:`0` at the beginning and end of the pulse.
 This is only meant to demonstrate the need for such an approach, and one simple example of one.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     from qiskit_dynamics import DiscreteSignal
     from qiskit_dynamics.signals import Convolution
@@ -134,7 +140,9 @@ This is only meant to demonstrate the need for such an approach, and one simple
 
 Observe, for example, the signal generated when all parameters are :math:`10^8`:
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     signal = signal_mapping(np.ones(80) * 1e8)
     signal.draw(t0=0., tf=signal.duration * signal.dt, n=1000, function='envelope')
@@ -148,7 +156,9 @@ the pulse via the standard fidelity measure:
 
 .. math:: f(U) = \frac{|\text{Tr}(XU)|^2}{4}
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     X_op = Operator.from_label('X').data
 
@@ -164,7 +174,9 @@ The function we want to optimize consists of:
 -  Simulating the Schrodinger equation over the length of the pulse sequence.
 -  Computing and return the infidelity (we minimize :math:`1 - f(U)`).
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     def objective(params):
 
@@ -199,7 +211,9 @@ Finally, we gradient optimize the objective:
 -  Call ``scipy.optimize.minimize`` with the above, with ``method='BFGS'`` and ``jac=True`` to
    indicate that the passed objective also computes the gradient.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     from jax import jit, value_and_grad
     from scipy.optimize import minimize
@@ -220,7 +234,9 @@ solver.
 We can draw the optimized signal, which is retrieved by applying the ``signal_mapping`` to the
 optimized parameters.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     opt_signal = signal_mapping(opt_results.x)
 
@@ -236,7 +252,9 @@ optimized parameters.
 Summing the signal samples yields approximately :math:`\pm 50`, which is equivalent to what one
 would expect based on a rotating wave approximation analysis.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     opt_signal.samples.sum()
 
@@ -252,7 +270,9 @@ instance, parameterized by ``sigma`` and ``width``. Although qiskit pulse provid
 :class:`~qiskit.pulse.library.GaussianSquare`, this class is not JAX compatible. See the user guide
 entry on :ref:`JAX-compatible pulse schedules <how-to use pulse schedules for jax-jit>`.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     import sympy as sym
     from qiskit import pulse
@@ -310,7 +330,9 @@ entry on :ref:`JAX-compatible pulse schedules <how-to use pulse schedules for ja
 Next, we construct a pulse schedule using the above parametrized Gaussian square pulse, convert it
 to a signal, and simulate the equation over the length of the pulse sequence.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     from qiskit_dynamics.pulse import InstructionToSignals
 
@@ -341,12 +363,16 @@ to a signal, and simulate the equation over the length of the pulse sequence.
 We set the initial values of ``sigma`` and ``width`` for the optimization as
 ``initial_params = np.array([10, 10])``.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     initial_params = np.array([10, 10])
     gaussian_square_generated_by_pulse(initial_params).draw()
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     from jax import jit, value_and_grad
     from scipy.optimize import minimize
@@ -367,6 +393,8 @@ We set the initial values of ``sigma`` and ``width`` for the optimization as
 
 We can draw the optimized pulse, whose parameters are retrieved by ``opt_results.x``.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     gaussian_square_generated_by_pulse(opt_results.x).draw()
\ No newline at end of file
diff --git a/docs/tutorials/qiskit_pulse.rst b/docs/tutorials/qiskit_pulse.rst
index 6b946265f..5d1d319d5 100644
--- a/docs/tutorials/qiskit_pulse.rst
+++ b/docs/tutorials/qiskit_pulse.rst
@@ -32,7 +32,9 @@ steps:
 
 First, we use the pulse module in Qiskit to create a pulse schedule.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     import numpy as np
     import qiskit.pulse as pulse
@@ -71,7 +73,9 @@ i.e. ``dt``, as well as the carrier frequency of the signals, i.e. ``w``. The pl
 envelopes and the signals resulting from this conversion. The dashed line shows the time at which
 the virtual ``Z`` gate is applied.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     from matplotlib import pyplot as plt
     from qiskit_dynamics.pulse import InstructionToSignals
@@ -99,7 +103,9 @@ to simulate pulse schedules. This requires specifying which channels act on whic
 carrier frequencies, and sample width ``dt``. Additionally, we setup this solver in the rotating
 frame and perform the rotating wave approximation.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     from qiskit.quantum_info.operators import Operator
     from qiskit_dynamics import Solver
@@ -131,17 +137,26 @@ In the last step we perform the simulation and plot the results. Note that, as w
 ``converter.get_signals`` above can also be passed to the ``signals`` argument and in this case
 should produce identical behavior.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
+    import time
     from qiskit.quantum_info.states import Statevector
 
     # Start the qubit in its ground state.
     y0 = Statevector([1., 0.])
 
-    %time sol = hamiltonian_solver.solve(t_span=[0., 2*T], y0=y0, signals=sxp, atol=1e-8, rtol=1e-8)
+    start_time = time.time()
+    
+    sol = hamiltonian_solver.solve(t_span=[0., 2*T], y0=y0, signals=sxp, atol=1e-8, rtol=1e-8)
 
+    print(f"Run time: {time.time() - start_time}")
 
-.. jupyter-execute::
+
+.. plot::
+    :context: close-figs
+    :include-source:
 
     def plot_populations(sol):
         pop0 = [psi.probabilities()[0] for psi in sol.y]
@@ -163,6 +178,8 @@ instruction in the pulse schedule. As expected, the first pulse moves the qubit
 the ``Y`` operator. Therefore, the second pulse, which drives around the ``Y``-axis due to the phase
 shift, has hardley any influence on the populations of the qubit.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     plot_populations(sol)
diff --git a/docs/userguide/how_to_configure_simulations.rst b/docs/userguide/how_to_configure_simulations.rst
index ac01486f0..697d57ac8 100644
--- a/docs/userguide/how_to_configure_simulations.rst
+++ b/docs/userguide/how_to_configure_simulations.rst
@@ -41,7 +41,9 @@ where:
 
 First, construct the components of the model:
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     import numpy as np
     from qiskit.quantum_info import Operator
@@ -70,7 +72,11 @@ First, construct the components of the model:
 Construct a :class:`.Solver` for the model as stated, without entering a rotating frame, and solve,
 timing the solver.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
+
+    import time
 
     solver = Solver(
         static_hamiltonian=static_hamiltonian,
@@ -78,12 +84,19 @@ timing the solver.
     )
 
     y0 = np.eye(dim, dtype=complex)
-    %time results = solver.solve(t_span=[0., T], y0=y0, signals=[drive_signal], atol=1e-10, rtol=1e-10)
+
+    start_time = time.time()
+
+    results = solver.solve(t_span=[0., T], y0=y0, signals=[drive_signal], atol=1e-10, rtol=1e-10)
+
+    print(f"Run time: {time.time() - start_time}")
 
 Next, define a :class:`.Solver` in the rotating frame of the static Hamiltonian by setting the
 ``rotating_frame`` kwarg, and solve, again timing the solver.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     rf_solver = Solver(
         static_hamiltonian=static_hamiltonian,
@@ -92,7 +105,12 @@ Next, define a :class:`.Solver` in the rotating frame of the static Hamiltonian
     )
 
     y0 = np.eye(dim, dtype=complex)
-    %time rf_results = rf_solver.solve(t_span=[0., T], y0=y0, signals=[drive_signal], atol=1e-10, rtol=1e-10)
+
+    start_time = time.time()
+
+    rf_results = rf_solver.solve(t_span=[0., T], y0=y0, signals=[drive_signal], atol=1e-10, rtol=1e-10)
+
+    print(f"Run time: {time.time() - start_time}")
 
 Observe that despite the two simulation problems being mathematically equivalent, it takes less time
 to solve in the rotating frame.
@@ -107,7 +125,9 @@ To compare the results, we use the fidelity function for unitary matrices:
 
 where :math:`d` is the dimension. A value of :math:`1` indicates equality of the unitaries.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     def fidelity(U, V):
         # the fidelity function
@@ -127,13 +147,17 @@ The discrepancy in solving times can be understood by examining the number of ri
 evaluations when solving the differential equation in each instance. The number of RHS evaluations
 for the first simulation (not in the rotating frame) was:
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     results.nfev
 
 Whereas the number of evaluations for the second simulation in the rotating frame was:
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     rf_results.nfev
 
@@ -152,7 +176,9 @@ RWA.
 Construct a solver for the same problem, now specifying an RWA cutoff frequency and the carrier
 frequencies relative to which the cutoff should be applied:
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     rwa_solver = Solver(
         static_hamiltonian=static_hamiltonian,
@@ -163,19 +189,28 @@ frequencies relative to which the cutoff should be applied:
     )
 
     y0 = np.eye(dim, dtype=complex)
-    %time rwa_results = rwa_solver.solve(t_span=[0., T], y0=y0, signals=[drive_signal], atol=1e-10, rtol=1e-10)
+
+    start_time = time.time()
+
+    rwa_results = rwa_solver.solve(t_span=[0., T], y0=y0, signals=[drive_signal], atol=1e-10, rtol=1e-10)
+
+    print(f"Run time: {time.time() - start_time}")
 
 We observe a further reduction in time, which is a result of the solver requiring even fewer RHS
 evaluations with the RWA:
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     rwa_results.nfev
 
 This speed comes at the cost of lower accuracy, owing to the fact that RWA is a legitimate
 *approximation*, which modifies the structure of the solution:
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     U_rwa = rwa_solver.model.rotating_frame.state_out_of_frame(T, rwa_results.y[-1])
 
@@ -205,7 +240,9 @@ Dynamics.
 
 Start off by configuring JAX.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     # configure jax to use 64 bit mode
     import jax
@@ -218,7 +255,9 @@ Reconstruct the model pieces at a much larger dimension, to observe the benefits
 arrays. Furthermore, set up the initial state to be a single column vector, to further highlight the
 benefits of the sparse representation.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     dim = 300
 
@@ -242,7 +281,9 @@ benefits of the sparse representation.
 Construct standard dense solver in the rotating frame of the static Hamiltonian, define a function
 to solve the system for a given amplitude, and just-in-time compile it using JAX.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     solver = Solver(
         static_hamiltonian=static_hamiltonian,
@@ -268,7 +309,9 @@ Construct sparse solver **in the frame of the diagonal of the static Hamiltonian
 function to solve the system for a given amplitude, and just-in-time compile it. Note that in this
 case the static Hamiltonian is already diagonal, but we explicitly highlight the need for this.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     sparse_solver = Solver(
         static_hamiltonian=static_hamiltonian,
@@ -293,21 +336,37 @@ case the static Hamiltonian is already diagonal, but we explicitly highlight the
 
 Run the dense simulation (twice to see the true compiled speed).
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
+
+    yf = jitted_dense_func(1.).block_until_ready()
+
+    start_time = time.time()
 
     yf = jitted_dense_func(1.).block_until_ready()
-    %time yf = jitted_dense_func(1.).block_until_ready()
+
+    print(f"Run time: {time.time() - start_time}")
 
 Run the sparse solver (twice to see the true compiled speed).
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
+
+    yf_sparse = jitted_sparse_func(1.).block_until_ready()
 
+    start_time = time.time()
+    
     yf_sparse = jitted_sparse_func(1.).block_until_ready()
-    %time yf_sparse = jitted_sparse_func(1.).block_until_ready()
+
+    print(f"Run time: {time.time() - start_time}")
 
 Verify equality of the results in a common frame.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     yf = solver.model.rotating_frame.state_out_of_frame(T, yf)
     yf_sparse = sparse_solver.model.rotating_frame.state_out_of_frame(T, yf_sparse)
diff --git a/docs/userguide/how_to_use_different_array_libraries.rst b/docs/userguide/how_to_use_different_array_libraries.rst
index e125f032b..df162ddfa 100644
--- a/docs/userguide/how_to_use_different_array_libraries.rst
+++ b/docs/userguide/how_to_use_different_array_libraries.rst
@@ -23,7 +23,9 @@ This guide addresses the following topics:
 
 First, configure JAX and import array libraries.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     # configure jax to use 64 bit mode
     import jax
@@ -38,7 +40,9 @@ First, configure JAX and import array libraries.
 
 Defining equivalent :class:`.Signal` instances, with envelope implemented in either NumPy or JAX.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     from qiskit_dynamics import Signal
 
@@ -54,19 +58,25 @@ Defining equivalent :class:`.Signal` instances, with envelope implemented in eit
 
 Evaluation of ``signal_numpy`` is executed with NumPy:
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     type(signal_numpy(0.1))
 
 Evaluation of ``signal_jax`` is executed with JAX:
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     type(signal_jax(0.1))
 
 JAX transformations can be applied to ``signal_jax``, e.g. just-in-time compilation:
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     from jax import jit
 
@@ -80,7 +90,9 @@ JAX transformations can be applied to ``signal_jax``, e.g. just-in-time compilat
 Internally, Qiskit Dynamics uses an extension of the default NumPy and SciPy array libraries offered
 by `Arraylias <https://qiskit-community.github.io/arraylias/>`_. These can be imported as:
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
     
     # alias for NumPy and corresponding aliased library
     from qiskit_dynamics import DYNAMICS_NUMPY_ALIAS
@@ -108,7 +120,9 @@ scans over a control parameter.
 
 First, we construct a :class:`.Solver` instance with a simple qubit model.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     import numpy as np
     from qiskit.quantum_info import Operator
@@ -136,7 +150,9 @@ Next, define the function to be compiled:
   - The output is the state of the system, starting in the ground state, at ``100`` points over the
     total evolution time.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     def sim_function(amp):
 
@@ -156,7 +172,9 @@ Next, define the function to be compiled:
 
 Compile the function.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     from jax import jit
     fast_sim = jit(sim_function)
@@ -165,22 +183,36 @@ The first time the function is called, JAX will compile an `XLA <https://www.ten
 version of the function, which is then executed. Hence, the time taken on the first call *includes*
 compilation time.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
-    %time ys = fast_sim(1.).block_until_ready()
+    start_time = time.time()
+
+    ys = fast_sim(1.).block_until_ready()
+
+    print(f"Run time: {time.time() - start_time}")
 
 
 On subsequent calls the compiled function is directly executed, demonstrating the true speed of the
 compiled function.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
+
+    start_time = time.time()
+    
+    fast_sim(1.).block_until_ready()
 
-    %timeit fast_sim(1.).block_until_ready()
+    print(f"Run time: {time.time() - start_time}")
 
 
 We use this function to plot the :math:`Z` expectation value over a range of input amplitudes.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     import matplotlib.pyplot as plt
 
diff --git a/docs/userguide/how_to_use_pulse_schedule_for_jax_jit.rst b/docs/userguide/how_to_use_pulse_schedule_for_jax_jit.rst
index 1c6223a44..0407cf627 100644
--- a/docs/userguide/how_to_use_pulse_schedule_for_jax_jit.rst
+++ b/docs/userguide/how_to_use_pulse_schedule_for_jax_jit.rst
@@ -30,7 +30,9 @@ Dynamics.
 
 First, configure JAX to run on CPU in 64 bit mode.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     # configure jax to use 64 bit mode
     import jax
@@ -48,7 +50,9 @@ Gaussian pulse to use in optimization, we need to instantiate a
 :class:`~qiskit.pulse.library.ScalableSymbolicPulse` with a Gaussian parameterization. First, define
 the symbolic representation in `sympy`.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     from qiskit import pulse
     from qiskit_dynamics.pulse import InstructionToSignals
@@ -75,7 +79,9 @@ the symbolic representation in `sympy`.
 
 Next, define the :class:`~qiskit.pulse.library.ScalableSymbolicPulse` using the above expression.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     _t, _duration, _amp, _sigma, _angle = sym.symbols("t, duration, amp, sigma, angle")
     _center = _duration / 2
@@ -105,7 +111,9 @@ Using a Gaussian pulse as an example, we show that a function involving
 :class:`~qiskit.pulse.library.ScalableSymbolicPulse` and the pulse to signal converter can be
 JAX-compiled (or more generally, JAX-transformed).
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     # use amplitude as the function argument
     def jit_func(amp):
diff --git a/docs/userguide/perturbative_solvers.rst b/docs/userguide/perturbative_solvers.rst
index 10d5ef489..941b7e921 100644
--- a/docs/userguide/perturbative_solvers.rst
+++ b/docs/userguide/perturbative_solvers.rst
@@ -53,7 +53,9 @@ First, configure JAX to run on CPU in 64 bit mode. See the :ref:`userguide on us
 different array libraries>` for a more detailed explanation of how to work with JAX in Qiskit
 Dynamics.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
     
     # configure jax to use 64 bit mode
     import jax
@@ -73,7 +75,9 @@ the model, which will slow down both the ODE solver and the initial construction
 solver. However after the initial construction, the higher frequencies in the model have no impact
 on the perturbative solver speed.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     import numpy as np
 
@@ -124,12 +128,15 @@ See the :class:`.DysonSolver` API docs for more details.
 
 For our example Hamiltonian we configure the :class:`.DysonSolver` as follows:
 
-.. jupyter-execute::
-
-    %%time
+.. plot::
+    :context: close-figs
+    :include-source:
 
+    import time
     from qiskit_dynamics import DysonSolver
 
+    start_time = time.time()
+
     dt = 0.1
     dyson_solver = DysonSolver(
         operators=[-1j * drive_hamiltonian],
@@ -143,6 +150,8 @@ For our example Hamiltonian we configure the :class:`.DysonSolver` as follows:
         rtol=1e-12
     )
 
+    print(f"Run time: {time.time() - start_time}")
+
 The above parameters are chosen so that the :class:`.DysonSolver` is fast and produces high accuracy
 solutions (measured and confirmed after the fact). The relatively large step size ``dt = 0.1`` is
 chosen for speed: the larger the step size, the fewer steps required. To ensure high accuracy given
@@ -159,7 +168,9 @@ over the interval ``[0, (T // dt) * dt]`` for an on-resonance drive with envelop
 ``envelope_func`` above. Running compiled versions of these functions gives a sense of the speeds
 attainable by these solvers.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     from qiskit_dynamics import Signal
     from jax import jit
@@ -180,17 +191,31 @@ attainable by these solvers.
 
 First run includes compile time.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
+
+    import time
+
+    start_time = time.time()
+    
+    yf_dyson = dyson_sim(1.).block_until_ready()
 
-    %time yf_dyson = dyson_sim(1.).block_until_ready()
+    print(f"Run time: {time.time() - start_time}")
 
 
 Once JIT compilation has been performance we can benchmark the performance of the JIT-compiled
 solver:
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
+
+    start_time = time.time()
+    
+    yf_dyson = dyson_sim(1.).block_until_ready()
 
-    %time yf_dyson = dyson_sim(1.).block_until_ready()
+    print(f"Run time: {time.time() - start_time}")
 
 
 4. Comparison to traditional ODE solver
@@ -199,7 +224,9 @@ solver:
 We now construct the same simulation using a standard solver to compare accuracy and simulation
 speed.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     from qiskit_dynamics import Solver
 
@@ -224,7 +251,9 @@ speed.
 
 Simulate with low tolerance for comparison to high accuracy solution.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     yf_low_tol = ode_sim(1., 1e-13)
     np.linalg.norm(yf_low_tol - yf_dyson)
@@ -232,22 +261,36 @@ Simulate with low tolerance for comparison to high accuracy solution.
 
 For speed comparison, compile at a tolerance with similar accuracy.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     jit_ode_sim = jit(lambda amp: ode_sim(amp, 1e-8))
 
-    %time yf_ode = jit_ode_sim(1.).block_until_ready()
+    start_time = time.time()
+
+    yf_ode = jit_ode_sim(1.).block_until_ready()
+
+    print(f"Run time: {time.time() - start_time}")
 
 Measure compiled time.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
-    %time yf_ode = jit_ode_sim(1.).block_until_ready()
+    start_time = time.time()
+
+    yf_ode = jit_ode_sim(1.).block_until_ready()
+
+    print(f"Run time: {time.time() - start_time}")
 
 
 Confirm similar accuracy solution.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     np.linalg.norm(yf_low_tol - yf_ode)
 
@@ -262,12 +305,14 @@ Next, we repeat our example using the Magnus-based perturbative solver. Setup of
 :class:`.MagnusSolver` is similar to the :class:`.DysonSolver`, but it uses the Magnus expansion and
 matrix exponentiation to simulate over each fixed time step.
 
-.. jupyter-execute::
-
-    %%time
+.. plot::
+    :context: close-figs
+    :include-source:
 
     from qiskit_dynamics import MagnusSolver
 
+    start_time = time.time()
+
     dt = 0.1
     magnus_solver = MagnusSolver(
         operators=[-1j * drive_hamiltonian],
@@ -281,10 +326,14 @@ matrix exponentiation to simulate over each fixed time step.
         rtol=1e-12
     )
 
+    print(f"Run time: {time.time() - start_time}")
+
 
 Setup simulation function.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     @jit
     def magnus_sim(amp):
@@ -299,18 +348,32 @@ Setup simulation function.
 
 First run includes compile time.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
-    %time yf_magnus = magnus_sim(1.).block_until_ready()
+    start_time = time.time()
+
+    yf_magnus = magnus_sim(1.).block_until_ready()
+
+    print(f"Run time: {time.time() - start_time}")
 
 Second run demonstrates speed of the simulation.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
+
+    start_time = time.time()
+    
+    yf_magnus = magnus_sim(1.).block_until_ready()
 
-    %time yf_magnus = magnus_sim(1.).block_until_ready()
+    print(f"Run time: {time.time() - start_time}")
 
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     np.linalg.norm(yf_magnus - yf_low_tol)
 
diff --git a/qiskit_dynamics/pulse/__init__.py b/qiskit_dynamics/pulse/__init__.py
index 99733ab58..d8a6382ea 100644
--- a/qiskit_dynamics/pulse/__init__.py
+++ b/qiskit_dynamics/pulse/__init__.py
@@ -40,8 +40,8 @@
 
 An example schedule, and the corresponding converted signals, is shown below.
 
-.. jupyter-execute::
-    :hide-code:
+.. plot::
+    :context: close-figs
 
     import matplotlib.pyplot as plt
     import matplotlib.gridspec as gridspec
diff --git a/qiskit_dynamics/signals/__init__.py b/qiskit_dynamics/signals/__init__.py
index 8e4e41225..9626ce515 100644
--- a/qiskit_dynamics/signals/__init__.py
+++ b/qiskit_dynamics/signals/__init__.py
@@ -142,8 +142,8 @@
 signal superimposed with sampled versions, both in the case of sampling the carrier, and in the case
 of sampling just the envelope (and keeping the carrier analog).
 
-.. jupyter-execute::
-    :hide-code:
+.. plot::
+    :context: close-figs
 
     from qiskit_dynamics.signals import Signal, DiscreteSignal
     from matplotlib import pyplot as plt
@@ -177,7 +177,9 @@
 carrier frequency of 400 MHz to create a signal at 500 MHz. Note that the code below does not make
 any assumptions about the time and frequency units which we interpret as ns and GHz, respectively.
 
-.. jupyter-execute::
+.. plot::
+    :context: close-figs
+    :include-source:
 
     import numpy as np
     from qiskit_dynamics.signals import DiscreteSignal, Sampler, IQMixer
diff --git a/requirements-dev.txt b/requirements-dev.txt
index 3516af256..86d98b2c5 100644
--- a/requirements-dev.txt
+++ b/requirements-dev.txt
@@ -2,7 +2,7 @@ stestr>=3.0.0
 astroid==2.14.2
 pylint==2.16.2
 black~=24.1
-qiskit-sphinx-theme~=1.16.0
+qiskit-sphinx-theme~=2.0
 sphinx-autodoc-typehints
 jupyter-sphinx
 # The following line is needed until