diff --git a/.ci/Project.toml b/.ci/Project.toml
index 5d94f0d8c..e1748090d 100644
--- a/.ci/Project.toml
+++ b/.ci/Project.toml
@@ -5,7 +5,6 @@ version = "0.1.0"
 
 [deps]
 Comonicon = "863f3e99-da2a-4334-8734-de3dacbe5542"
-CondaPkg = "992eb4ea-22a4-4c89-a5bb-47a3300528ab"
 CoverageTools = "c36e975a-824b-4404-a568-ef97ca766997"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 LiveServer = "16fef848-5104-11e9-1b77-fb7a48bbb589"
diff --git a/.ci/src/example.jl b/.ci/src/example.jl
index bbd60f5e8..006ab6898 100644
--- a/.ci/src/example.jl
+++ b/.ci/src/example.jl
@@ -141,7 +141,6 @@ in parallel.
     example_dir = root_dir("examples")
     script = """
     using Pkg
-    using CondaPkg
     using Literate
     for name in readdir(\"$example_dir\")
         project_dir = joinpath(\"$example_dir\", name)
@@ -149,7 +148,6 @@ in parallel.
 
         Pkg.activate(project_dir)
         Pkg.instantiate()
-        CondaPkg.resolve()
 
         @info "building" project_dir
         Literate.$target(
diff --git a/CondaPkg.toml b/CondaPkg.toml
deleted file mode 100644
index 255866961..000000000
--- a/CondaPkg.toml
+++ /dev/null
@@ -1,2 +0,0 @@
-[deps]
-matplotlib = "3.5.1"
diff --git a/Project.toml b/Project.toml
index 4beb2b54d..428673598 100644
--- a/Project.toml
+++ b/Project.toml
@@ -9,11 +9,11 @@ BloqadeLattices = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe4"
 BloqadeMIS = "bd27d05e-4ce1-5e74-84dd-c5d7d508bbe2"
 BloqadeODE = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe5"
 BloqadeWaveforms = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe7"
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7"
-PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Yao = "5872b779-8223-5990-8dd0-5abbb0748c8c"
 YaoSubspaceArrayReg = "bd27d05e-4ce1-5e79-84dd-c5d7d508ade2"
@@ -29,7 +29,6 @@ ColorSchemes = "3"
 Colors = "0.12"
 ForwardDiff = "0.10"
 Measurements = "2"
-PythonCall = "0.8,0.9"
 Reexport = "1"
 Yao = "0.9"
 YaoSubspaceArrayReg = "0.2"
diff --git a/docs/Project.toml b/docs/Project.toml
index a879c23f1..2e879db8c 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -16,7 +16,6 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
-PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 Yao = "5872b779-8223-5990-8dd0-5abbb0748c8c"
 YaoAPI = "0843a435-28de-4971-9e8b-a9641b2983a8"
diff --git a/docs/src/install.md b/docs/src/install.md
index 10fa5fb18..2231bd2c0 100644
--- a/docs/src/install.md
+++ b/docs/src/install.md
@@ -57,8 +57,6 @@ Built on top of Ubuntu Server 20.04 LTS, this image includes
 - [Yao.jl](https://yaoquantum.org/)
 - [Revise.jl](https://github.com/timholy/Revise.jl)
 - [BenchmarkTools.jl](https://juliaci.github.io/BenchmarkTools.jl/stable/)
-- [PythonCall.jl](https://cjdoris.github.io/PythonCall.jl/stable/) 
-- Conda package manager, provided by [Miniconda](https://docs.conda.io/en/latest/miniconda.html) 
 
 #### Bloqade CUDA AMI
 
@@ -227,8 +225,6 @@ Built on top of Ubuntu Server 20.04 LTS, this image includes
 - [Yao.jl](https://yaoquantum.org/)
 - [Revise.jl](https://github.com/timholy/Revise.jl)
 - [BenchmarkTools.jl](https://juliaci.github.io/BenchmarkTools.jl/stable/)
-- [PythonCall.jl](https://cjdoris.github.io/PythonCall.jl/stable/)
-- Conda package manager, provided by [Mamba](https://mamba.readthedocs.io/en/latest/index.html)
 - Jupyter Lab interface with dedicated Julia and Python kernels
 - Integrated Terminal for interactive command-line sessions
 
diff --git a/docs/src/waveform.md b/docs/src/waveform.md
index 7d3134fbd..c02c8f8db 100644
--- a/docs/src/waveform.md
+++ b/docs/src/waveform.md
@@ -17,14 +17,11 @@ which is a created by providing a callable object and a real number `duration`:
 Bloqade gives users the flexibility to specify general waveforms by inputting functions. The following code constructs a sinusoidal waveform with a time duration of 2 μs:
 
 ```@example waveform
-using Bloqade
-using PythonCall
-plt = pyimport("matplotlib.pyplot")
+using Bloqade, Bloqade.CairoMakie
 waveform = Waveform(t->2.2*2π*sin(2π*t), duration = 2);
 Bloqade.plot(waveform)
 ```
-In our documentation, we use the
-python package [`matplotlib`](https://matplotlib.org) for plotting.
+In our documentation, we use the [`CairoMakie`](https://docs.makie.org/) for plotting.
 
 Bloqade supports built-in waveforms for convenience (see References below). 
 For example, the codes below create different waveform shapes with a single line:
@@ -63,7 +60,9 @@ wf1 = piecewise_linear(;clocks, values=values1);
 values2 = 2π*rand(length(clocks)-1)
 wf2 = piecewise_constant(;clocks, values=values2); 
 
-fig, (ax1, ax2) = plt.subplots(figsize=(12, 4), ncols=2)
+fig = Figure(size=(960, 400))
+ax1 = Axis(fig[1, 1])
+ax2 = Axis(fig[1, 2])
 Bloqade.plot!(ax1, wf1)
 Bloqade.plot!(ax2, wf2)
 fig
@@ -104,7 +103,9 @@ waveform:
 wf = piecewise_linear(clocks=[0.0, 2.0, 3.0, 4.0], values=2π*[0.0, 3.0, 1.1, 2.2]);
 swf = smooth(wf;kernel_radius=0.1);
 
-fig, (ax1, ax2) = plt.subplots(figsize=(12, 4), ncols=2)
+fig = Figure(size=(960, 400))
+ax1 = Axis(fig[1, 1])
+ax2 = Axis(fig[1, 2])
 Bloqade.plot!(ax1, wf)
 Bloqade.plot!(ax2, swf)
 fig
@@ -120,7 +121,9 @@ wf2 = sinusoidal(duration = 2.2, amplitude = 2.2*2π);
 wf3 = wf1 + wf2; 
 wf4 = wf1 - wf2;
 
-fig, (ax1, ax2) = plt.subplots(figsize=(12, 4), ncols=2)
+fig = Figure(size=(960, 400))
+ax1 = Axis(fig[1, 1])
+ax2 = Axis(fig[1, 2])
 Bloqade.plot!(ax1, wf3)
 Bloqade.plot!(ax2, wf4)
 fig
@@ -133,7 +136,9 @@ To increase the strength of a waveform by some factors, we can directly use `*`:
 wf = linear_ramp(;duration=2.2, start_value=0.0, stop_value=1.0*2π);
 wf_t = 3 * wf;
 
-fig, (ax1, ax2) = plt.subplots(figsize=(12, 4), ncols=2)
+fig = Figure(size=(960, 400))
+ax1 = Axis(fig[1, 1])
+ax2 = Axis(fig[1, 2])
 Bloqade.plot!(ax1, wf)
 Bloqade.plot!(ax2, wf_t)
 fig
@@ -144,7 +149,9 @@ Such operations can also be broadcasted by using `.*`:
 ```@example waveform
 wf2, wf3 = [2.0, 3.0] .* wf1; 
 
-fig, (ax1, ax2) = plt.subplots(figsize=(12, 4), ncols=2)
+fig = Figure(size=(960, 400))
+ax1 = Axis(fig[1, 1])
+ax2 = Axis(fig[1, 2])
 Bloqade.plot!(ax1, wf2)
 Bloqade.plot!(ax2, wf3)
 fig
diff --git a/examples/1.blockade/Project.toml b/examples/1.blockade/Project.toml
index 6de2daf32..d0d02f70a 100644
--- a/examples/1.blockade/Project.toml
+++ b/examples/1.blockade/Project.toml
@@ -7,7 +7,6 @@ BloqadeMIS = "bd27d05e-4ce1-5e74-84dd-c5d7d508bbe2"
 BloqadeODE = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe5"
 BloqadeWaveforms = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe7"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
-PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 Yao = "5872b779-8223-5990-8dd0-5abbb0748c8c"
 YaoArrayRegister = "e600142f-9330-5003-8abb-0ebd767abc51"
diff --git a/examples/1.blockade/main.jl b/examples/1.blockade/main.jl
index 251b84bc5..17079e1ec 100644
--- a/examples/1.blockade/main.jl
+++ b/examples/1.blockade/main.jl
@@ -345,26 +345,29 @@ for _ in TimeChoiceIterator(integrator2, 0.0:dt:Tmax)
 end
 
 # Plot the data:
-using PythonCall # Use matplotlib to generate plots
-matplotlib = pyimport("matplotlib")
-plt = pyimport("matplotlib.pyplot")
-
-ax = plt.subplot(1, 1, 1)
-plt.plot(times, real(densities), "k", label = "Full space")
-plt.plot(times, real(densities2), "r--", label = "Subspace")
-ax.axis([0, Tmax, 0, 0.45])
-plt.xlabel("Time (us)")
-plt.ylabel("Rydberg density")
-plt.tight_layout()
-plt.legend()
-
-inset_axes = pyimport("mpl_toolkits.axes_grid1.inset_locator")
-ax2 = inset_axes.inset_axes(ax, width = "20%", height = "30%", loc = "lower right", borderpad = 1)
-plt.plot(times, real(densities - densities2))
-plt.axis([0, 0.5, -0.001, 0.003])
-plt.ylabel("Difference", fontsize = 12)
-plt.yticks(LinRange(-0.001, 0.003, 5), fontsize = 12);
-plt.xticks([0, 0.2, 0.4, 0.6], fontsize = 12);
+using Bloqade.CairoMakie # Use CairoMakie to generate plots
+
+fig = Figure()
+ax = Axis(fig[1, 1], xlabel = "Time (us)", ylabel = "Rydberg density", limits=(0.0, Tmax, 0.0, 0.45))
+lines!(ax, times, real(densities), label = "Full space", color=:black)
+lines!(ax, times, real(densities2), label = "Subspace", color=:red, linestyle=:dash)
+axislegend(ax, position = :rt)
+fig
+
+# The inset axis
+ax2 = Axis(fig[1, 1],
+    width=Relative(0.2),
+    height=Relative(0.3),
+    halign=0.9,
+    valign=0.1,
+    limits=(0.0, 0.5, -0.001, 0.003),
+    ylabel="Difference",
+    yticks=LinRange(-0.001, 0.003, 5),
+    xticks=[0, 0.2, 0.4, 0.6],
+    backgroundcolor=:lightgray)
+
+lines!(ax2, times, real(densities - densities2), color = :black)
+fig;
 
 # ![RydbergBlockadeSubspace](../../../assets/RydbergBlockadeSubspace.png)
 
diff --git a/examples/2.adiabatic/Project.toml b/examples/2.adiabatic/Project.toml
index 604de30b8..6a814e65b 100644
--- a/examples/2.adiabatic/Project.toml
+++ b/examples/2.adiabatic/Project.toml
@@ -7,7 +7,6 @@ BloqadeMIS = "bd27d05e-4ce1-5e74-84dd-c5d7d508bbe2"
 BloqadeODE = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe5"
 BloqadeWaveforms = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe7"
 KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
-PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d"
 YaoSubspaceArrayReg = "bd27d05e-4ce1-5e79-84dd-c5d7d508ade2"
 
 [extras]
diff --git a/examples/2.adiabatic/main.jl b/examples/2.adiabatic/main.jl
index 109502f71..bfaab0c4a 100644
--- a/examples/2.adiabatic/main.jl
+++ b/examples/2.adiabatic/main.jl
@@ -19,12 +19,10 @@
 # Let's start by importing the required libraries:
 
 using Bloqade
-using PythonCall
+using Bloqade.CairoMakie
 using KrylovKit
 using SparseArrays
 
-plt = pyimport("matplotlib.pyplot");
-
 # # Ground State Properties
 
 # We start by probing the ground state properties of the Rydberg Hamiltonian in a 1D system. 
@@ -59,21 +57,15 @@ end
 
 # To compare, we first plot the density profile when ``\Delta= -2π \times 10`` MHz: 
 
-fig, ax = plt.subplots(figsize = (10, 4))
-ax.bar(1:nsites, density_g[1, :])
-ax.set_xticks(1:nsites)
-ax.set_xlabel("Sites")
-ax.set_ylabel("Rydberg density")
-ax.set_title("Density Profile: 1D Chain, Δ = -2π * 10 MHz")
+fig = Figure(size=(800, 320))
+ax = Axis(fig[1, 1], xlabel = "Sites", ylabel = "Rydberg density", title = "Density Profile: 1D Chain, Δ = -2π * 10 MHz", xticks = 1:nsites)
+barplot!(ax, 1:nsites, density_g[1, :])
 fig
 
 # We can see that the Rydberg densities in this case are close to 0 for all sites. In contrast, for ``\Delta= 2π \times 10`` MHz, the density shows a clear ``Z_2`` ordered profile:
-fig, ax = plt.subplots(figsize = (10, 4))
-ax.bar(1:nsites, density_g[30, :])
-ax.set_xticks(1:nsites)
-ax.set_xlabel("Sites")
-ax.set_ylabel("Rydberg density")
-ax.set_title("Density Profile: 1D Chain, Δ = 2π * 10 MHz")
+fig = Figure(size=(800, 320))
+ax = Axis(fig[1, 1], xlabel = "Sites", ylabel = "Rydberg density", title = "Density Profile: 1D Chain, Δ = 2π * 10 MHz", xticks = 1:nsites)
+barplot!(ax, 1:nsites, density_g[30, :])
 fig
 
 # More generally, we can plot an order parameter as a function of ``\Delta`` to clearly see the onset of phase transition. 
@@ -83,13 +75,10 @@ order_para = map(1:Δ_step) do ii
     return sum(density_g[ii, 1:2:nsites]) - sum(density_g[ii, 2:2:nsites]) # density on odd sites - density on even sites
 end
 
-fig, ax = plt.subplots(figsize = (10, 4))
-ax.plot(Δ / 2π, order_para)
-ax.set_xlabel("Δ/2π (MHz) ")
-ax.set_ylabel("Order parameter")
+fig = Figure(size=(800, 320))
+ax = Axis(fig[1, 1], xlabel = "Δ/2π (MHz)", ylabel = "Order parameter")
+lines!(ax, Δ / (2π), order_para)
 fig
-
-# From the density profile of ground states and the change in the order parameter, we can observe a phase transition with changing ``\Delta``. 
 # Below, we show that by slowly changing the parameters of the Hamiltonian, we can follow the trajectory of the ground states and adiabatically evolve the atoms from the ground state to the ``Z_2`` 
 # ordered state.
 
@@ -106,13 +95,12 @@ U2 = 2π * 10;
 Δ = piecewise_linear(clocks = [0.0, 0.6, 2.1, total_time], values = [U1, U1, U2, U2]);
 
 # We plot the two waveforms:
-fig, (ax1, ax2) = plt.subplots(ncols = 2, figsize = (12, 4))
-Bloqade.plot!(ax1, Ω)
-ax1.set_ylabel("Ω/2π (MHz)")
-Bloqade.plot!(ax2, Δ)
-ax2.set_ylabel("Δ/2π (MHz)")
+fig = Figure(size=(960, 320))
+ax1 = Axis(fig[1, 1], ylabel = "Ω/2π (MHz)")
+ax2 = Axis(fig[1, 2], ylabel = "Δ/2π (MHz)")
+Bloqade.plot!(ax1, Ω / (2π))
+Bloqade.plot!(ax2, Δ / (2π))
 fig
-
 # We generate the positions for a 1D atomic chain again: 
 
 nsites = 9
@@ -136,19 +124,13 @@ densities = []
 for _ in TimeChoiceIterator(integrator, 0.0:1e-3:total_time)
     push!(densities, rydberg_density(reg))
 end
-D = hcat(densities...);
-
-# and finally plot the time-dependent dynamics of Rydberg density for each site:
-fig, ax = plt.subplots(figsize = (10, 4))
-shw = ax.imshow(real(D), interpolation = "nearest", aspect = "auto", extent = [0, total_time, 0.5, nsites + 0.5])
-ax.set_xlabel("time (μs)")
-ax.set_ylabel("site")
-ax.set_xticks(0:0.2:total_time)
-ax.set_yticks(1:nsites)
-bar = fig.colorbar(shw)
-fig
+D = hcat(densities...)';
 
-# We can clearly see that a ``Z_2`` ordered state has been generated by the specified adiabatic pulse sequence. 
+fig = Figure(size=(800, 320))
+ax = Axis(fig[1, 1], xlabel = "time (μs)", ylabel = "site", xticks=0:0.2:total_time, yticks=(0.5:1:nsites, string.(1:nsites)))
+shw = heatmap!(ax, LinRange(0.0, total_time, size(D, 1)), LinRange(0.5, nsites - 0.5, size(D, 2)), real(D))
+Colorbar(fig[1, 2], shw, label = "Rydberg density")
+fig
 # We can also confirm it by plotting the bitstring distribution at the final time step:
 
 bitstring_hist(reg; nlargest = 20)
@@ -194,11 +176,11 @@ total_time = 2.9
 U = 2π * 15.0
 Δ = piecewise_linear(clocks = [0.0, 0.3, 2.6, total_time], values = [-U, -U, U, U]);
 
-fig, (ax1, ax2) = plt.subplots(ncols = 2, figsize = (10, 4))
+fig = Figure(size=(800, 320))
+ax1 = Axis(fig[1, 1], ylabel = "Ω/2π (MHz)")
+ax2 = Axis(fig[1, 2], ylabel = "Δ/2π (MHz)")
 Bloqade.plot!(ax1, Ω)
-ax1.set_ylabel("Ω/2π (MHz)")
 Bloqade.plot!(ax2, Δ)
-ax2.set_ylabel("Δ/2π (MHz)")
 fig
 
 # Then, we use the waveforms and atom positions to create a Hamiltonian and define a time evolution problem:
@@ -213,13 +195,10 @@ densities = [];
 for _ in TimeChoiceIterator(integrator, 0.0:1e-3:total_time)
     push!(densities, rydberg_density(reg))
 end
-D = hcat(densities...)
-
-fig, ax = plt.subplots(figsize = (10, 4))
-shw = ax.imshow(real(D), interpolation = "nearest", aspect = "auto", extent = [0, total_time, 0.5, nsites + 0.5])
-ax.set_xlabel("time (μs)")
-ax.set_ylabel("site")
-ax.set_xticks(0:0.2:total_time)
-ax.set_yticks(1:nsites)
-bar = fig.colorbar(shw)
+D = hcat(densities...)'
+
+fig = Figure(size=(800, 320))
+ax = Axis(fig[1, 1], xlabel = "time (μs)", ylabel = "site", xticks=0:0.2:total_time, yticks=(0.5:1:nsites, string.(1:nsites)))
+shw = heatmap!(ax, LinRange(0.0, total_time, size(D, 1)), LinRange(0.5, nsites - 0.5, size(D, 2)), real(D))
+Colorbar(fig[1, 2], shw, label = "Rydberg density")
 fig
diff --git a/examples/3.quantum-scar/Project.toml b/examples/3.quantum-scar/Project.toml
index 9c8f6f348..6b49bc0fa 100644
--- a/examples/3.quantum-scar/Project.toml
+++ b/examples/3.quantum-scar/Project.toml
@@ -1,5 +1,4 @@
 [deps]
-CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 Bloqade = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe1"
 BloqadeExpr = "bd27d05e-4ce1-5e79-84dd-c5d7d508abe2"
 BloqadeKrylov = "bd27d05e-4cd1-5e79-84dd-c5d7d508ade2"
@@ -7,8 +6,8 @@ BloqadeLattices = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe4"
 BloqadeMIS = "bd27d05e-4ce1-5e74-84dd-c5d7d508bbe2"
 BloqadeODE = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe5"
 BloqadeWaveforms = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe7"
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
-PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 Yao = "5872b779-8223-5990-8dd0-5abbb0748c8c"
 YaoArrayRegister = "e600142f-9330-5003-8abb-0ebd767abc51"
diff --git a/examples/3.quantum-scar/main.jl b/examples/3.quantum-scar/main.jl
index f172a7871..bbb0922df 100644
--- a/examples/3.quantum-scar/main.jl
+++ b/examples/3.quantum-scar/main.jl
@@ -20,10 +20,9 @@
 # To start, we first import the required libraries:
 
 using Bloqade
-using PythonCall
+using Bloqade.CairoMakie
 using Random
 
-plt = pyimport("matplotlib.pyplot");
 # # Many-Body Rabi Oscillations with Rydberg Blockade
 
 # We first demonstrate that the strong Rydberg interactions have important effects on the Rabi oscillations of Rydberg atoms.
@@ -57,11 +56,9 @@ end
 # Here, we use a [`KrylovEvolution`](@ref) object to simulate the dynamics for a time-independent Hamiltonian.
 # One can also use ODE to simulate the dynamics. For an example, see the [Adiabatic Evolution](@ref) tutorial.
 # The Rydberg density of this atom exhibits Rabi oscillations as a function of time, shown by the plot below:
-fig, ax = plt.subplots()
-ax.plot(clocks, density1[1, :])
-ax.set_xlabel("Time (μs)")
-ax.set_ylabel("Single Rydberg Probability")
-ax.set_title("Rabi Oscillation: Single Atom Case")
+fig = Figure()
+ax = Axis(fig[1,1], xlabel = "Time (μs)", ylabel = "Single Rydberg Probability", title = "Rabi Oscillation: Single Atom Case")
+lines!(ax, clocks, density1[1, :], label = "1 atom")
 fig
 
 # For the case of 2 and 3 atoms, if they are separated far enough with negligible interactions, the total Rydberg excitation densities are simply the sum of the Rydberg density for each atom.
@@ -92,14 +89,12 @@ density3 = sum(density3, dims = 1);
 # where ``N`` is the number of atoms.
 # For more information, please refer to [H. Bernien, et al. (10.1038/nature24622)](https://www.nature.com/articles/nature24622).
 # The total Rydberg density for the 1-, 2-, and 3-atom system is plotted below:
-fig, ax = plt.subplots()
-ax.plot(clocks, density1[1, :])
-ax.plot(clocks, density2[1, :])
-ax.plot(clocks, density3[1, :])
-ax.set_xlabel("Time (μs)")
-ax.set_ylabel("Rydberg Probability")
-ax.set_title("Many-body Rabi Oscillation for 1-, 2-, and 3-atom system")
-ax.legend(["1 atom", "2 atoms", "3 atoms"], loc = "lower right")
+fig = Figure()
+ax = Axis(fig[1,1], xlabel = "Time (μs)", ylabel = "Rydberg Probability", title = "Many-body Rabi Oscillation for 1-, 2-, and 3-atom system")
+lines!(ax, clocks, density1[1, :], label = "1 atom")
+lines!(ax, clocks, density2[1, :], label = "2 atoms")
+lines!(ax, clocks, density3[1, :], label = "3 atoms")
+axislegend(position = :rb)
 fig
 
 # From this plot, we can see that the total Rydberg density for 2 (3) atom case does not exceed 1. This is because
@@ -134,11 +129,13 @@ atoms = generate_sites(ChainLattice(), nsites, scale = 5.72)
 Ω_tot = append(Ω1, Ω2);
 Δ_tot = append(Δ1, Δ2);
 
-fig, (ax1, ax2) = plt.subplots(ncols = 2, figsize = (12, 4))
+fig = Figure(size=(960, 320))
+ax1 = Axis(fig[1, 1], ylabel = "Ω/2π (MHz)")
+ax2 = Axis(fig[1, 2], ylabel = "Δ/2π (MHz)")
 Bloqade.plot!(ax1, Ω_tot)
 Bloqade.plot!(ax2, Δ_tot)
-ax1.set_ylabel("Ω/2π (MHz)")
-ax2.set_ylabel("Δ/2π (MHz)")
+fig[1, 1] = ax1
+fig[1, 2] = ax2
 fig
 
 # Note that the total evolution time is 4.2 μs.
@@ -171,28 +168,23 @@ end
 # We first plot the Rydberg density for each site as a function of time:
 
 clocks = 0:1e-3:total_time
-D = hcat(densities...)
-
-fig, ax = plt.subplots(figsize = (10, 4))
-shw = ax.imshow(real(D), interpolation = "nearest", aspect = "auto", extent = [0, total_time, 0.5, nsites + 0.5])
-ax.set_xlabel("time (μs)")
-ax.set_ylabel("site")
-ax.set_xticks(0:0.4:total_time)
-ax.set_yticks(1:nsites)
-bar = fig.colorbar(shw)
-fig
+D = hcat(densities...)'
 
+fig = Figure(size=(800, 320))
+ax = Axis(fig[1, 1]; xlabel = "Time (μs)", ylabel = "Site", xticks = 0:0.4:total_time, yticks = (0.5:1:nsites, string.(1:nsites)))
+shw = heatmap!(ax, LinRange(0, total_time, size(D, 1)), LinRange(0.5, nsites - 0.5, size(D, 2)), real(D), colormap = :viridis)
+bar = Colorbar(fig[1, 2], shw)
+fig
 # We can see that the state evolves to a Neel state after the first part of the pulse (time around 2.2 μs). 
 # After that, there are clear oscillations between the two patterns of the Rydberg density, which is a signature of the quantum scar.
 
 # We can also plot the entanglement as a function of time:
 
-fig, ax = plt.subplots(figsize = (10, 4))
-ax.plot(clocks, entropy)
-ax.set_xlabel("time (μs)")
-ax.set_ylabel("entanglement entropy")
+fig = Figure(size=(800, 320))
+ax = Axis(fig[1, 1]; xlabel = "time (μs)", ylabel = "entanglement entropy")
+line = lines!(ax, clocks, entropy)
+fig[1, 1] = ax
 fig
-
 # # A Different Initial State 
 
 # In order to see that the revivals depends strongly on the initial state, 
@@ -204,26 +196,18 @@ clocks = 0.0:1e-2:total_time;
 
 init_d = product_state(bit"100000101")
 prob_d = KrylovEvolution(init_d, clocks, hd)
-density_mat_d = zeros(nsites, length(clocks))
+density_mat_d = zeros(length(clocks), nsites)
 
 for info in prob_d
     for i in 1:nsites
-        density_mat_d[i, info.step] = rydberg_density(info.reg, i)
+        density_mat_d[info.step, i] = rydberg_density(info.reg, i)
     end
 end
 
-fig, ax = plt.subplots(figsize = (10, 4))
-shw = ax.imshow(
-    real(density_mat_d),
-    interpolation = "nearest",
-    aspect = "auto",
-    extent = [0, total_time, 0.5, nsites + 0.5],
-)
-ax.set_xlabel("time (μs)")
-ax.set_ylabel("site")
-ax.set_xticks(0:0.2:total_time)
-ax.set_yticks(1:nsites)
-bar = fig.colorbar(shw)
+fig = Figure(size=(800, 320))
+ax = Axis(fig[1, 1]; xlabel = "time (μs)", ylabel = "site", xticks = 0:0.2:total_time, yticks = (0.5:1.0:nsites, string.(1:nsites)))
+shw = heatmap!(ax, LinRange(0.0, total_time, size(density_mat_d, 1)), LinRange(0.5, nsites - 0.5, size(density_mat_d, 2)), density_mat_d, colormap = :viridis)
+bar = Colorbar(fig[1, 2], shw)
 fig
 
 # From this figure, we see that the density does not show long-lived oscillations.
diff --git a/examples/4.LGT/Project.toml b/examples/4.LGT/Project.toml
index 604de30b8..6a814e65b 100644
--- a/examples/4.LGT/Project.toml
+++ b/examples/4.LGT/Project.toml
@@ -7,7 +7,6 @@ BloqadeMIS = "bd27d05e-4ce1-5e74-84dd-c5d7d508bbe2"
 BloqadeODE = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe5"
 BloqadeWaveforms = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe7"
 KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
-PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d"
 YaoSubspaceArrayReg = "bd27d05e-4ce1-5e79-84dd-c5d7d508ade2"
 
 [extras]
diff --git a/examples/4.LGT/main.jl b/examples/4.LGT/main.jl
index 227fbdac5..e9fca0173 100644
--- a/examples/4.LGT/main.jl
+++ b/examples/4.LGT/main.jl
@@ -13,9 +13,7 @@
 # We first import the required packages 
 
 using Bloqade
-using PythonCall
-plt = pyimport("matplotlib.pyplot");
-
+using Bloqade.CairoMakie
 
 # ## Mapping between the Rydberg system and the LGT
 
@@ -59,16 +57,13 @@ total_time = 3.5;
 Ω1 = piecewise_linear(clocks = [0.0, 0.2, total_time], values = [0.0, Ωmax, Ωmax]);
 
 # The waveforms for the detuning and the Rabi frequency are shown below
-fig1, (ax1, ax2) = plt.subplots(ncols = 2, figsize = (12, 4))
+fig1 = Figure(size=(960, 320))
+ax1 = Axis(fig1[1, 1], ylabel="Ω1/2π (MHz)", xgridvisible=true, ygridvisible=true)
+ax2 = Axis(fig1[1, 2], ylabel="Δ1/2π (MHz)", xgridvisible=true, ygridvisible=true)
 Bloqade.plot!(ax1, Ω1)
 Bloqade.plot!(ax2, Δ1)
-ax1.set_ylabel("Ω1/2π (MHz)")
-ax2.set_ylabel("Δ1/2π (MHz)")
-ax1.grid()
-ax2.grid()
 fig1
 
-
 # In order to simulate the gauge theory dynamics, we define the function
 # `get_average_rydberg_densities` which takes in a detuning
 # and Rabi frequency $(Δ, Ω)$, and returns the final state and the Rydberg density: 
@@ -93,11 +88,11 @@ end;
 # by simulating the dynamics governed by the waveforms, followed by plotting the density profile, as shown below:
 
 dens1 = get_average_rydberg_densities(Δ1, Ω1) ;
-fig2, ax = plt.subplots(figsize = (10, 4)) ;
-ax.bar(1:N, dens1[end]) ;
-ax.set_xlabel("site index")
-ax.set_ylabel("Average Rydberg densities")
-fig2 
+
+fig2 = Figure(size = (800, 320))
+ax = Axis(fig2[1, 1], xlabel="site index", ylabel="Average Rydberg densities")
+barplot!(ax, 1:N, dens1[end])
+fig2
 
 # Recalling the mapping between the Rydberg chain and the LGT illustrated above, 
 # we see that the final state is an approximation of the anti-string state of the LGT, or a $Z_2$ ordered state of the Rydberg chain. 
@@ -148,7 +143,10 @@ end ;
 
 # As an example, for the case with a single defect, we show the detuning for the central site, which is the defect, and those for other sites separately below:
 
-fig3, (ax1, ax2) = plt.subplots(ncols = 2, figsize = (12, 4))
+fig3 = Figure(size=(960, 320))
+ax1 = Axis(fig3[1, 1], title="Detuning for the central site", xgridvisible=true, ygridvisible=true)
+ax2 = Axis(fig3[1, 2], title="Detunings for other sites", xgridvisible=true, ygridvisible=true)
+
 for idx in 1 : length(atoms)
     if idx == floor(Int, N/2)+1
         Bloqade.plot!(ax1, Δ2_single_defect[idx])
@@ -156,10 +154,6 @@ for idx in 1 : length(atoms)
         Bloqade.plot!(ax2, Δ2_single_defect[idx])
     end
 end
-ax1.grid()
-ax2.grid()
-ax1.set_title("Detuning for the central site")
-ax2.set_title("Detunings for other sites")
 
 fig3
 
@@ -169,14 +163,13 @@ fig3
 dens2 = get_average_rydberg_densities(Δ2_single_defect, Ω2)
 dens3 = get_average_rydberg_densities(Δ2_two_defects, Ω2)
 
-fig4, (ax1, ax2) = plt.subplots(nrows = 2, figsize = (10, 4), frameon=false)
-ax1.bar(1:N, dens2[end])
-ax2.bar(1:N, dens3[end])
-fig4.supxlabel("site index")
-fig4.supylabel("Average Rydberg densities", x=0.06)
+fig4 = Figure(size=(800, 320))
+ax1 = Axis(fig4[1, 1], xlabel="site index", ylabel="Average Rydberg densities")
+ax2 = Axis(fig4[2, 1], xlabel="site index", ylabel="Average Rydberg densities")
+barplot!(ax1, 1:N, dens2[end])
+barplot!(ax2, 1:N, dens3[end])
 fig4
 
-
 # Again, we see that the Rydberg density at the defects are not exactly zero, but the prepared states,
 # as we shall see below, serve as good initial states to study the propagation of particle-antiparticle pairs in LGT. 
 
@@ -197,20 +190,15 @@ end
 # Again, as an example, for the case with a single defect, we show the detuning for the central site, 
 # which is the defect, and those for other sites separately below:
 
-fig5, (ax1, ax2) = plt.subplots(ncols = 2, figsize = (12, 4))
-for idx in 1 : length(atoms)
-    if idx == floor(Int, N/2)+1
-        Bloqade.plot!(ax1, Δ3_single_defect[idx])
-    else
-        Bloqade.plot!(ax2, Δ3_single_defect[idx])
-    end
-end
-Bloqade.plot!(ax1, Ω3)
-ax1.grid()
-ax2.grid()
-ax1.legend(["Detuning for the central site", "Rabi frequency for all sites"])
-ax2.legend(["Detunings for other sites"], loc="center left")
+fig5 = Figure(size=(960, 320))
+ax1 = Axis(fig5[1, 1], xlabel="Time", ylabel="Detuning", xgridvisible=true, ygridvisible=true)
+ax2 = Axis(fig5[1, 2], xlabel="Time", ylabel="Detuning", xgridvisible=true, ygridvisible=true)
 
+Bloqade.plot!(ax1, Δ3_single_defect[floor(Int, N/2)+1], label="Detuning for the central site")
+Bloqade.plot!(ax2, Δ3_single_defect[1], label="Detunings for other sites")
+Bloqade.plot!(ax1, Ω3; label="Rabi frequency for all sites")
+axislegend(ax1, position=:rb)
+axislegend(ax2, position=:rb)
 fig5
 
 # ### Simulating Particle-Antiparticle Pairs in LGT Dynamics 
@@ -235,16 +223,15 @@ clocks = clocks[ind0: end];
 # Then we plot the Rydberg density as a function of time, where the two panels correspond to the cases with single and two defects respectively:
 yticks = range(clocks[1], stop=clocks[end], length=10);
 yticks = [string(ytick)[1:4] for ytick in yticks][end:-1:1];
-fig6, (ax1, ax2) = plt.subplots(ncols = 2, figsize = (12, 8), sharey=true)
-ax1.imshow(transpose(D_single_defect)[end:-1:1,:], aspect="auto", interpolation="none")
-ax1.set_xlabel("sites")
-ax1.set_ylabel("time (μs)")
-ax1.set_yticks(range(1, stop = length(clocks), length=10), yticks)
-im = ax2.imshow(transpose(D_two_defects)[end:-1:1,:], aspect="auto", interpolation="none")
-ax2.set_xlabel("sites")
-fig6.colorbar(im, ax=[ax1, ax2])
-fig6 
 
+fig6 = Figure(size=(960, 640))
+ax1 = Axis(fig6[1, 1], xlabel = "sites", ylabel = "time (μs)", yticks = (range(1, stop = length(clocks), length = 10), string.(round.(range(clocks[1], stop=clocks[end], length=10); digits=2))))
+ax2 = Axis(fig6[1, 2], xlabel = "sites", yticklabelsvisible=false)
+clims = (0.0, 1.0)
+img1 = image!(ax1, D_single_defect[end:-1:1, :], colormap = :inferno, colorrange=clims, interpolate=false)
+img2 = image!(ax2, D_two_defects[end:-1:1, :], colormap = :inferno, colorrange=clims, interpolate=false)
+Colorbar(fig6[1, 3]; limits=clims)
+fig6
 # From the left panel, we can observe a light-cone-shaped region originating from the particle-antiparticle pair in the vacuum. 
 # At the right panel, we show the interference of two light cones, which produces an additional change of periodicity corresponding to the elastic scattering. 
 # When the particle or antiparticle reaches the boundary of the chain, it will be scattered back as observed. 
diff --git a/examples/5.MIS/Project.toml b/examples/5.MIS/Project.toml
index c345c5027..16ee9abaa 100644
--- a/examples/5.MIS/Project.toml
+++ b/examples/5.MIS/Project.toml
@@ -10,7 +10,6 @@ BloqadeWaveforms = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe7"
 GenericTensorNetworks = "3521c873-ad32-4bb4-b63d-f4f178f42b49"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
-PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 YaoSubspaceArrayReg = "bd27d05e-4ce1-5e79-84dd-c5d7d508ade2"
 
diff --git a/examples/5.MIS/main.jl b/examples/5.MIS/main.jl
index 418a74042..8266bb856 100644
--- a/examples/5.MIS/main.jl
+++ b/examples/5.MIS/main.jl
@@ -34,8 +34,7 @@ using Bloqade
 using Random
 using GenericTensorNetworks
 using Optim
-using PythonCall
-plt = pyimport("matplotlib.pyplot");
+using Bloqade.CairoMakie
 
 # ## Setting Up the Problem
 
@@ -59,7 +58,8 @@ Bloqade.plot(atoms, blockade_radius = 7.5)
 # the [generic tensor network algorithm (J. Liu et al. (10.48550/arXiv.2205.03718))](https://arxiv.org/pdf/2205.03718.pdf) 
 # in the package [`GenericTensorNetworks`](https://github.com/QuEraComputing/GenericTensorNetworks.jl):
 graph = BloqadeMIS.unit_disk_graph(atoms, 7.5)
-mis_size_and_counting = GenericTensorNetworks.solve(IndependentSet(graph), CountingMax())[]
+independent_set_problem = IndependentSet(graph)
+mis_size_and_counting = GenericTensorNetworks.solve(GenericTensorNetwork(independent_set_problem), CountingMax())[]
 
 # The `solve` function takes a graph instance and a solution space property as inputs,
 # where the graph instance is generated by the [`unit_disk_graph`](@ref) function in the Bloqade submodule `BloqadeMIS`,
@@ -80,11 +80,11 @@ T_max = 0.6
 Δ_end = 2π * 11
 Δ = piecewise_linear(clocks = [0.0, 0.1, 0.5, T_max], values = [Δ_start, Δ_start, Δ_end, Δ_end])
 
-fig, (ax1, ax2) = plt.subplots(ncols = 2, figsize = (12, 4))
+fig = Figure(size=(960, 320))
+ax1 = Axis(fig[1, 1]; ylabel="Ω/2π (MHz)")
+ax2 = Axis(fig[1, 2]; ylabel="Δ/2π (MHz)")
 Bloqade.plot!(ax1, Ω)
-ax1.set_ylabel("Ω/2π (MHz)")
 Bloqade.plot!(ax2, Δ)
-ax2.set_ylabel("Δ/2π (MHz)")
 fig
 
 # Here, the total time is fixed to `T_max`, the adiabatic evolution path is specified by the [`piecewise_linear`](@ref) function.
@@ -111,7 +111,7 @@ bitstring_hist(prob.reg; nlargest = 20)
 # The correctness of the output can be verified by comparing it to the classical solution:
 
 best_bit_strings = most_probable(prob.reg, 2)
-all_optimal_configs = GenericTensorNetworks.solve(IndependentSet(graph), ConfigsMax())[]
+all_optimal_configs = GenericTensorNetworks.solve(GenericTensorNetwork(independent_set_problem), ConfigsMax())[]
 @assert all(bs -> GenericTensorNetworks.StaticBitVector([bs...]) ∈ all_optimal_configs.c, best_bit_strings)
 
 # We can also visualize these atoms and check them visually:
@@ -175,11 +175,12 @@ clocks = [0, cumsum(durations)...]
 Ω2 = piecewise_constant(; clocks = clocks, values = repeat([Ω_max, 0.0], 3))
 Δ2 = piecewise_constant(; clocks = clocks, values = repeat([0.0, Δ_end], 3))
 
-fig, (ax1, ax2) = plt.subplots(ncols = 2, figsize = (12, 4))
+
+fig = Figure(size=(960, 320))
+ax1 = Axis(fig[1, 1]; ylabel="Ω/2π (MHz)")
+ax2 = Axis(fig[1, 2]; ylabel="Δ/2π (MHz)")
 Bloqade.plot!(ax1, Ω2)
-ax1.set_ylabel("Ω/2π (MHz)")
 Bloqade.plot!(ax2, Δ2)
-ax2.set_ylabel("Δ/2π (MHz)")
 fig
 
 # The `piecewise_constant` pulses can be more accurately simulated with the [`KrylovEvolution`](@ref) solver.
@@ -231,8 +232,8 @@ end
 
 # Let us check the loss function again using the initial point above: 
 x0 = durations
-rydberg_density, reg1 = loss_piecewise_constant(atoms, x0)
-rydberg_density
+density, reg1 = loss_piecewise_constant(atoms, x0)
+density
 
 # The most probable bitstrings are:
 
@@ -267,11 +268,11 @@ bitstring_hist(reg1_final; nlargest = 20)
 pulse_piecewise_linear = piecewise_linear(clocks = [0.0, 0.05, 0.1, 0.5, 0.55, T_max], values = [0, 0, 0.4, 0.4, 0, 0]);
 pulse_smooth = smooth(pulse_piecewise_linear; kernel_radius = 0.02);
 
-fig, ax = plt.subplots()
-Bloqade.plot!(ax, pulse_piecewise_linear)
-Bloqade.plot!(ax, pulse_smooth)
-ax.set_ylabel("strength")
-ax.legend(["piecewise linear", "smoothened piecewise linear"], loc = "lower right")
+fig = Figure()
+ax = Axis(fig[1, 1]; ylabel="strength")
+Bloqade.plot!(ax, pulse_piecewise_linear; label = "piecewise linear")
+Bloqade.plot!(ax, pulse_smooth; label = "smoothened piecewise linear")
+axislegend(ax; position=:rb)
 fig
 
 # Here, the function [`smooth`](@ref) takes a `kernel_radius` keyword parameter as the Gaussian kernel parameter.
@@ -317,15 +318,15 @@ T_max = 0.6
     values = [Δ_start, Δ_start, Δ0 * x0[1], Δ0 * x0[2], Δ0 * x0[3], Δ_end, Δ_end],
 )
 
-rydberg_density, reg2, Δ_initial_smooth = loss_piecewise_linear(atoms, x0)
-rydberg_density
+density, reg2, Δ_initial_smooth = loss_piecewise_linear(atoms, x0)
+density
 
 # and plot the smoothened waveform:
-fig, ax = plt.subplots()
-Bloqade.plot!(ax, Δ_initial)
-Bloqade.plot!(ax, Δ_initial_smooth)
-ax.set_ylabel("Δ/2π (MHz)")
-ax.legend(["piecewise linear", "smoothened piecewise linear"], loc = "lower right")
+fig = Figure()
+ax = Axis(fig[1, 1]; ylabel="Δ/2π (MHz)")
+Bloqade.plot!(ax, Δ_initial; label = "piecewise linear")
+Bloqade.plot!(ax, Δ_initial_smooth; label = "smoothened piecewise linear")
+axislegend(ax; position=:rb)
 fig
 
 # Let's plot the distribution:
@@ -344,11 +345,11 @@ bitstring_hist(reg_final; nlargest = 20)
 
 # We can also plot out the final optimized waveform for Δ 
 # and compare it with the initial waveform:
-fig, ax = plt.subplots()
-Bloqade.plot!(ax, Δ_initial_smooth)
-Bloqade.plot!(ax, Δ_final)
-ax.set_ylabel("Δ/2π (MHz)")
-ax.legend(["initial", "optimized"], loc = "lower right")
+fig = Figure()
+ax = Axis(fig[1, 1]; ylabel="Δ/2π (MHz)")
+Bloqade.plot!(ax, Δ_initial_smooth; label = "initial")
+Bloqade.plot!(ax, Δ_final; label = "optimized")
+axislegend(ax; position=:rb)
 fig
 
 # ## Applications in VLSI Chip Manufacturing 
@@ -380,35 +381,28 @@ fig
 using Bloqade
 using GenericTensorNetworks
 using GenericTensorNetworks: unit_disk_graph
-using PythonCall
+using Bloqade.CairoMakie
 using Random
 using Optim
-plt = pyimport("matplotlib.pyplot");
-patches = pyimport("matplotlib.patches");
-np = pyimport("numpy")
 
 # ### Setting up the Problem
 
 # We can visualize this problem as a grid of 64K RAM chips, where some of the grid cells are removed. 
 
 
-fig, ax = plt.subplots(figsize=(5, 5))
-
 G = 5 # grid size (G x G)
 defects = []
 for i in range(1, G*G//15)
     push!(defects, ((rand(Int)%G+10) % G, (rand(Int)%G + 10) % G))
 end
 
-ax.grid(visible=true)
-plt.xticks(1:G) 
-plt.yticks(1:G) 
-
+fig = Figure(size=(400, 400))
+ax = Axis(fig[1, 1]; ygridvisible=true, xgridvisible=true, xticks=1:G, yticks=1:G, aspect=DataAspect(), limits=(0, G, 0, G))
 grid = zeros((G, G))
 
 function plot_defects(defects)
     for (x, y) in defects
-        ax.add_patch(patches.Rectangle((x, y), 1, 1, color="red"))
+        CairoMakie.poly!(Rect(x, y, 1.0, 1.0), color = :red)
         grid[x+1, y+1] = 1
     end
 end
@@ -446,9 +440,9 @@ for x in 1:G-1
 end
 
 R_opt = sqrt(2sqrt(2))/2
+circles = []
 for circ in circs
-    patch = patches.Circle(circ, R_opt, fill=false)
-    ax.add_patch(patch)
+    push!(circles, CairoMakie.poly!(Circle(CairoMakie.Point(Float64.(circ)), R_opt), color = :transparent, strokewidth=1))
 end
 
 fig
@@ -476,7 +470,8 @@ Bloqade.plot(atoms, blockade_radius = R)
 # We will first use classical algorithms to find the optimal solution to the problem:
 
 unit_graph = unit_disk_graph(atoms, R)
-configs = GenericTensorNetworks.solve(IndependentSet(unit_graph), ConfigsMax())[]
+independent_set_problem = IndependentSet(unit_graph)
+configs = GenericTensorNetworks.solve(GenericTensorNetwork(independent_set_problem), ConfigsMax())[]
 MIS_config = configs.c[1]
 
 # We can visualize the solution on the atom lattice:
@@ -485,24 +480,10 @@ Bloqade.plot(atoms, blockade_radius = R; colors = [iszero(b) ? "white" : "red" f
 
 # We can also visualize the solution on the chip grid:
 
-function clear_circ()
-    for i in 1:5
-        for patch in ax.patches
-            try
-                patch.radius
-                patch.remove()
-            catch
-                print(patch)
-            end
-        end
-    end
-end
-
-clear_circ()
+delete!.(Ref(ax.scene), circles)
 for (circ, used) in zip(circs, MIS_config)
     if used == 1
-        patch = patches.Circle(circ, R_opt, fill=false)
-        ax.add_patch(patch)
+        CairoMakie.poly!(Circle(CairoMakie.Point(Float64.(circ)), R_opt), color = :transparent, strokewidth=1)
     end
 end
 
@@ -521,12 +502,12 @@ T_max = 0.6
 Δ_end = 2π * 11
 Δ = piecewise_linear(clocks = [0.0, 0.1, 0.5, T_max], values = [Δ_start, Δ_start, Δ_end, Δ_end])
 
-plot, (ax1, ax2) = plt.subplots(ncols = 2, figsize = (12, 4))
+fig = Figure(size=(960, 320))
+ax1 = Axis(fig[1, 1]; ylabel="Ω/2π (MHz)")
+ax2 = Axis(fig[1, 2]; ylabel="Δ/2π (MHz)")
 Bloqade.plot!(ax1, Ω)
-ax1.set_ylabel("Ω/2π (MHz)")
 Bloqade.plot!(ax2, Δ)
-ax2.set_ylabel("Δ/2π (MHz)")
-plot
+fig
 
 # And construct the Hamiltonian:
 
@@ -586,12 +567,12 @@ end
 function plot_waves(params::Vector{Float64})
     Ω, Δ = get_waves(params)
 
-    graph, (ax1, ax2) = plt.subplots(ncols = 2, figsize = (12, 4))
+    fig = Figure(size=(960, 320))
+    ax1 = Axis(fig[1, 1]; ylabel="Ω/2π (MHz)")
+    ax2 = Axis(fig[1, 2]; ylabel="Δ/2π (MHz)")
     Bloqade.plot!(ax1, Ω)
-    ax1.set_ylabel("Ω/2π (MHz)")
     Bloqade.plot!(ax2, Δ)
-    ax2.set_ylabel("Δ/2π (MHz)")
-    return graph
+    return fig
 end
 
 x0 = [0.1, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8]
diff --git a/examples/6.MWIS/Project.toml b/examples/6.MWIS/Project.toml
index a20f21a7a..34714e734 100644
--- a/examples/6.MWIS/Project.toml
+++ b/examples/6.MWIS/Project.toml
@@ -11,6 +11,5 @@ GenericTensorNetworks = "3521c873-ad32-4bb4-b63d-f4f178f42b49"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
-PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 YaoSubspaceArrayReg = "bd27d05e-4ce1-5e79-84dd-c5d7d508ade2"
diff --git a/examples/6.MWIS/main.jl b/examples/6.MWIS/main.jl
index 5f498fe01..d3199258f 100644
--- a/examples/6.MWIS/main.jl
+++ b/examples/6.MWIS/main.jl
@@ -16,9 +16,7 @@ Random.seed!(42)
 using Graphs
 using GenericTensorNetworks
 using GenericTensorNetworks: unit_disk_graph
-using Bloqade
-using PythonCall
-plt = pyimport("matplotlib.pyplot");
+using Bloqade, Bloqade.CairoMakie
 
 # We now specify the atom locations and construct an example  
 # unit disk graph on a square lattice with nearest-neighbor connections. 
@@ -26,7 +24,7 @@ plt = pyimport("matplotlib.pyplot");
 # and all vertices closer than a distance 1.5 are connected by an edge:
 locs = [(1, -1), (4, 0), (1, 1), (2, 0), (0, 0), (2, 2), (2, -2), (3, 1), (3, -1)];
 g = unit_disk_graph(locs, 1.5)
-show_graph(g; locs = locs, vertex_colors = ["white" for i in 1:nv(g)])
+show_graph(g, map(x->60 .* x, locs); vertex_colors = ["white" for i in 1:nv(g)])
 
 # We then assign random weights to each vertex for this example problem:
 weights = [rand() for i in 1:nv(g)];
@@ -34,9 +32,10 @@ weights = [rand() for i in 1:nv(g)];
 # We can solve the MWIS problem classically for this graph
 # using the [GenericTensorNetworks](https://github.com/QuEraComputing/GenericTensorNetworks.jl) package. 
 # The MWIS is shown in red.
-configs_mapped = GenericTensorNetworks.solve(IndependentSet(g; weights = weights), ConfigsMax())[]
+wmis_problem = IndependentSet(g, weights)
+configs_mapped = GenericTensorNetworks.solve(GenericTensorNetwork(wmis_problem), ConfigsMax())[]
 MIS_config = configs_mapped.c[1]
-show_graph(g; locs = locs, vertex_colors = [iszero(MIS_config[i]) ? "white" : "red" for i in 1:nv(g)])
+show_graph(g, map(x->60 .* x, locs), vertex_colors = [iszero(MIS_config[i]) ? "white" : "red" for i in 1:nv(g)])
 
 # # Quantum Adiabatic Algorithm to Solve the MWIS Problem
 
@@ -85,13 +84,13 @@ end
 t_max = 1.5
 H, Ω, Δ = build_adiabatic_sweep(g, Ω_max, Δ_max, t_max, weights);
 
-fig, (ax1, ax2) = plt.subplots(ncols = 2, figsize = (14, 4))
+fig = Figure(; size=(700, 200))
+ax1 = Axis(fig[1, 1]; ylabel = "Ω/2π (MHz)")
+ax2 = Axis(fig[1, 2]; ylabel = "Δ/2π (MHz)")
 Bloqade.plot!(ax1, Ω)
-ax1.set_ylabel("Ω/2π (MHz)")
 for i in 1:nv(g)
     Bloqade.plot!(ax2, Δ[i])
 end
-ax2.set_ylabel("Δ/2π (MHz)")
 fig
 
 # ## Compute the MIS Probability and the Adiabatic Timescale
@@ -130,14 +129,12 @@ a, b = linear_fit(t_list[P_MWIS.>0.9], y)
 T_LZ = -1 / b;
 
 # Finally, we plot the results:
-fig, (ax1, ax2) = plt.subplots(ncols = 2, figsize = (16, 6))
-ax1.scatter(t_list, P_MWIS)
-ax1.set_ylabel("MWIS Probability")
-ax1.set_xlabel("Time (μs)")
-
-ax2.scatter(t_list, broadcast(log, 1 .- P_MWIS))
-ax2.plot(t_list, a .+ b .* t_list)
-ax2.set_xlabel("Time (μs)")
-ax2.set_ylabel("log(1 - MWIS Probability)")
+fig = Figure(; size=(800, 300))
+ax1 = Axis(fig[1, 1]; ylabel = "MWIS Probability", xlabel = "Time (μs)")
+ax2 = Axis(fig[1, 2]; ylabel = "log(1 - MWIS Probability)", xlabel = "Time (μs)")
+
+scatter!(ax1, t_list, P_MWIS)
+scatter!(ax2, t_list, broadcast(log, 1 .- P_MWIS))
+lines!(ax2, t_list, a .+ b .* t_list)
 
 fig
diff --git a/examples/7.QMC/Project.toml b/examples/7.QMC/Project.toml
index 199fb4b6d..3c5182df2 100644
--- a/examples/7.QMC/Project.toml
+++ b/examples/7.QMC/Project.toml
@@ -1,10 +1,17 @@
 [deps]
 BinningAnalysis = "b7192094-8e58-5052-a244-180a858778ee"
 Bloqade = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe1"
+BloqadeExpr = "bd27d05e-4ce1-5e79-84dd-c5d7d508abe2"
+BloqadeKrylov = "bd27d05e-4cd1-5e79-84dd-c5d7d508ade2"
+BloqadeLattices = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe4"
+BloqadeMIS = "bd27d05e-4ce1-5e74-84dd-c5d7d508bbe2"
+BloqadeODE = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe5"
 BloqadeQMC = "2c0eea08-e05e-11ec-33df-b7e33af02477"
+BloqadeWaveforms = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe7"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Yao = "5872b779-8223-5990-8dd0-5abbb0748c8c"
+YaoSubspaceArrayReg = "bd27d05e-4ce1-5e79-84dd-c5d7d508ade2"
diff --git a/examples/8.simulating noise/Project.toml b/examples/8.simulating noise/Project.toml
index 9a0ca79f7..93bc05154 100644
--- a/examples/8.simulating noise/Project.toml	
+++ b/examples/8.simulating noise/Project.toml	
@@ -1,6 +1,12 @@
 [deps]
 Bloqade = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe1"
+BloqadeExpr = "bd27d05e-4ce1-5e79-84dd-c5d7d508abe2"
+BloqadeKrylov = "bd27d05e-4cd1-5e79-84dd-c5d7d508ade2"
+BloqadeLattices = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe4"
+BloqadeMIS = "bd27d05e-4ce1-5e74-84dd-c5d7d508bbe2"
 BloqadeNoisy = "7534646f-0cda-4b9e-a311-3a9166d831c9"
+BloqadeODE = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe5"
+BloqadeWaveforms = "bd27d05e-4ce1-5e79-84dd-c5d7d508bbe7"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
@@ -11,3 +17,4 @@ Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 ProgressBars = "49802e3a-d2f1-5c88-81d8-b72133a6f568"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 Yao = "5872b779-8223-5990-8dd0-5abbb0748c8c"
+YaoSubspaceArrayReg = "bd27d05e-4ce1-5e79-84dd-c5d7d508ade2"
diff --git a/src/Bloqade.jl b/src/Bloqade.jl
index 8509fadd1..da92fb7f9 100644
--- a/src/Bloqade.jl
+++ b/src/Bloqade.jl
@@ -42,15 +42,7 @@ using Reexport
 
 export rydberg_density, rydberg_corr, bitstring_hist, bitstring_hist!, get_average_rydberg_densities
 
-using PythonCall
-const plt = PythonCall.pynew()
-
-function __init__()
-    # copied from PyPlotCall.jl
-    PythonCall.pycopy!(plt, pyimport("matplotlib.pyplot"))
-    return
-end
-
+using CairoMakie
 using Colors, ColorSchemes
 
 include("plots/plots.jl")
diff --git a/src/plots/bitstring_hist.jl b/src/plots/bitstring_hist.jl
index 30f387756..3878c5993 100644
--- a/src/plots/bitstring_hist.jl
+++ b/src/plots/bitstring_hist.jl
@@ -11,27 +11,24 @@ end
 """
     bitstring_hist!(ax, register; nlargest::Int, title="", kw...)
 
-Plot the bitstring histgram.
+Plot the bitstring histogram.
 
 # Arguments
 
-- `ax`: the axis object from `matplotlib.pyplot`.
+- `ax`: the axis object from `CairoMakie`.
 - `register`: the register to plot.
 
 # Keyword Arguments
 
 - `nlargest`: plot the first `nlargest` bitstrings.
 - `title`: title of the plot.
-- `kw`: other keyword supported by `matplotlib.bar` function.
+- `kw`: other keyword arguments supported by `CairoMakie.bars!` function.
 """
 function bitstring_hist!(ax, r::ArrayReg; nlargest::Int, title = "", kw...)
     ps = probs(r)
     indices = find_largest(ps, nlargest)
 
-    obj = ax.bar(1:length(indices), ps[indices], kw...)
-    ax.set_xlabel("bitstring")
-    ax.set_ylabel("probability")
-    ax.set_xticks(1:length(indices), string.(indices .- 1; base = 2, pad = nqubits(r)), rotation = 60, ha = "right")
+    obj = barplot!(ax, 1:length(indices), ps[indices]; kw...)
     return obj
 end
 
@@ -40,22 +37,14 @@ function bitstring_hist!(ax, r::SubspaceArrayReg; nlargest::Int, title = "", kw.
     probs = abs2.(state)
     indices = find_largest(probs, nlargest)
 
-    obj = ax.bar(1:length(indices), probs[indices]; kw...)
-    ax.set_xlabel("bitstring")
-    ax.set_ylabel("probability")
-    ax.set_xticks(
-        1:length(indices),
-        string.(space(r).subspace_v[indices]; base = 2, pad = nqubits(r)),
-        rotation = 60,
-        ha = "right",
-    )
+    obj = barplot!(ax, 1:length(indices), probs[indices]; kw...)
     return obj
 end
 
 """
     bitstring_hist(r; kw...)
 
-Plot the bitstring histgram.
+Plot the bitstring histogram.
 
 # Arguments
 
@@ -65,10 +54,24 @@ Plot the bitstring histgram.
 
 - `nlargest`: plot the first `nlargest` bitstrings.
 - `title`: title of the plot.
-- `kw...`: other keyword supported by `matplotlib.bar` function.
+- `kw...`: other keyword arguments supported by `CairoMakie.bars!` function.
 """
-function bitstring_hist(r; kw...)
-    fig, ax = plt.subplots()
-    bitstring_hist!(ax, r; kw...)
+function bitstring_hist(r; nlargest::Int, kw...)
+    fig = CairoMakie.Figure()
+    xticks_string = _xticks_string(r, nlargest)
+    ax = Axis(fig[1, 1], xlabel="bitstring", ylabel="probability", xticks = (1:length(xticks_string), xticks_string), xticklabelrotation=π/3)
+
+    bitstring_hist!(ax, r; nlargest, kw...)
     return fig
 end
+function _xticks_string(r::ArrayReg, nlargest::Int)
+    ps = probs(r)
+    indices = find_largest(ps, nlargest)
+    return string.(indices .- 1; base = 2, pad = nqubits(r))
+end
+function _xticks_string(r::SubspaceArrayReg, nlargest::Int)
+    state = statevec(r)
+    probs = abs2.(state)
+    indices = find_largest(probs, nlargest)
+    return string.(space(r).subspace_v[indices]; base = 2, pad = nqubits(r))
+end
\ No newline at end of file
diff --git a/src/plots/waveform.jl b/src/plots/waveform.jl
index 4ec37a309..d305dea9b 100644
--- a/src/plots/waveform.jl
+++ b/src/plots/waveform.jl
@@ -1,13 +1,12 @@
-function plot!(ax, wf::Waveform)
+function plot!(ax, wf::Waveform; kwargs...)
     clocks = sample_clock(wf)
-    fig = ax.plot(clocks, BloqadeWaveforms._rm_err.(sample_values(wf, clocks) ./ (2π));)
-    ax.set_xlabel("time (μs)")
-    ax.set_ylabel("value / 2π (MHz)")
+    lines!(ax, clocks, BloqadeWaveforms._rm_err.(sample_values(wf, clocks) ./ (2π)); kwargs...)
     return ax
 end
 
-function plot(wf::Waveform)
-    fig, ax = plt.subplots()
-    plot!(ax, wf)
+function plot(wf::Waveform; kwargs...)
+    fig = CairoMakie.Figure()
+    ax = Axis(fig[1, 1]; xlabel="time (μs)", ylabel = "value / 2π (MHz)")
+    plot!(ax, wf; kwargs...)
     return fig
-end
+end
\ No newline at end of file