diff --git a/docs/Project.toml b/docs/Project.toml
index 5f8dc1f..9ca1c02 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,5 +1,6 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+DemoCards = "311a05b2-6137-4a5a-b473-18580a3d38b5"
 PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
 Batsrus = "e74ebddf-6ac1-4047-a0e5-c32c99e57753"
 
diff --git a/docs/examples/basics/config.json b/docs/examples/basics/config.json
new file mode 100644
index 0000000..3b4c5c6
--- /dev/null
+++ b/docs/examples/basics/config.json
@@ -0,0 +1,7 @@
+{
+   "theme": "nocoverlist",
+   "order": [
+     "demo_animate_1d.md",
+     "demo_animate_2d.md"
+   ]
+}
\ No newline at end of file
diff --git a/docs/examples/basics/demo_animate_1d.md b/docs/examples/basics/demo_animate_1d.md
new file mode 100644
index 0000000..6f021af
--- /dev/null
+++ b/docs/examples/basics/demo_animate_1d.md
@@ -0,0 +1,106 @@
+# ---
+# title: Line Animation
+# id: demo_line_animation
+# date: 2024-05-08
+# author: "[Hongyang Zhou](https://github.com/henry2004y)"
+# julia: 1.10.3
+# description: 1D line animation using pyplot
+# ---
+
+This example shows how to create 1D space line animation from series of SWMF outputs.
+
+* Vertical line moving at the solar wind speed for reference.
+* Twin axes sharing the same x-axis for showing two quantities in one frame.
+* Time-dependent title.
+
+```julia
+using Batsrus, PyPlot, Printf
+
+"""
+    animate1d(files::Vector{String}, vars::Vector{String}; kwargs...)
+
+Saving series of plots of `vars` from SWMF output `files`.
+
+# Keywords
+- `filedir::String="./"`: input SWMF file directory.
+- `outdir::String="out/"`: output directory.
+- `output_interval::Int=1`: Timestep interval for output files.
+- `vmin=-1`: plot value lower bound.
+- `vmax=1`: plot value upper bound.
+- `refloc=1e5`: reference vertical line initial location.
+"""
+function animate1d(files::Vector{String}, vars::Vector{String};
+   filedir::String="./", outdir::String="out/", output_interval::Int=1,
+   vmin=-1, vmax=1, refloc=1e5)
+   nfile = length(files)
+   x = let
+      bd = load(joinpath(filedir, files[1]))
+      bd.x[:,1,1]
+   end
+
+   fig = plt.figure(figsize=(12,5), constrained_layout=true)
+   ax = plt.axes(xlim=extrema(x), ylim=(vmin, vmax))
+
+   line, = ax.plot([], [], lw=1)
+   # Move reference line forward in the canvas
+   vl = ax.axvline(0, ls="-", color="tab:brown", lw=1, zorder=10)
+
+   color = "tab:blue"
+   ax.set_xlabel("x [km]"; fontsize=14)
+   ax.set_ylabel(var*" [nT]"; fontsize=14, color)
+   ax.tick_params(axis="y", labelcolor=color)
+
+   ax2 = ax.twinx()
+   ax2.set_ylim(-600, 0)
+
+   color = "tab:red"
+   ax2.set_ylabel("Ux [km/s]"; fontsize=14, color)
+
+   line2, = ax2.plot([], []; color, alpha=0.7)
+   ax2.tick_params(axis="y", labelcolor=color)
+
+
+   for (i, file) in enumerate(files)
+      @info "$i in $nfile"
+      outname = outdir*lpad(i, 4, '0')*".png"
+      isfile(outname) && continue
+
+      bd = load(joinpath(filedir, file))
+
+      d = bd[vars[1]][:,1]
+      title_str = @sprintf "t = %4.1f s" bd.head.time
+      line.set_data(x, d)
+      d2 = bd[vars[2]][:,1]
+      line2.set_data(x, d2)
+
+      ax.set_title(title_str)
+
+      refloc -= output_interval * VSW
+      if refloc <= 0
+         refloc += 1e5
+      end
+      vl.set_xdata([refloc, refloc])
+
+      savefig(outname, bbox_inches="tight", dpi=200)
+   end
+
+   close()
+end
+
+####################
+# Constants
+const VSW = 500.0  # Solar wind speed, [km/s]
+
+# Data directory
+filedir = "./"
+# Plot variables
+vars = ["By", "uxs1"]
+
+# Find simulation data
+files = let
+   pick_file = file -> startswith(file, "z") && endswith(file, ".out")
+   filter(pick_file, readdir(filedir))
+end
+
+animate1d(files, vars)
+```
diff --git a/docs/examples/basics/demo_animate_2d.md b/docs/examples/basics/demo_animate_2d.md
new file mode 100644
index 0000000..abbc15e
--- /dev/null
+++ b/docs/examples/basics/demo_animate_2d.md
@@ -0,0 +1,151 @@
+# ---
+# title: Contour Animation
+# id: demo_contour_animation
+# date: 2024-05-08
+# author: "[Hongyang Zhou](https://github.com/henry2004y)"
+# julia: 1.10.3
+# description: 2D contour animation using pyplot
+# ---
+
+This example shows how to create 2D colored contour animation from series of SWMF outputs.
+
+* Time-dependent title.
+* Tweakable colorbar.
+* Plotting streamline on top of colored contours is a bit hacky because there is currently no intrinsic methods for removing the streamlines. Our solution here is to dispatch a streamline into lines and arrows and then remove them separately.
+
+```julia
+using Batsrus, PyPlot, Printf
+
+"""
+    animate2d(files::Vector{String}, var::String; kwargs...)
+
+Saving series of plots of `var` from SWMF output `files`.
+
+# Keywords
+- `filedir::String="./"`: input SWMF file directory.
+- `outdir::String="out/"`: output directory.
+- `vmin=-4`: plot value lower bound.
+- `vmax=4`: plot value upper bound.
+- `cmap=matplotlib.cm.RdBu_r`: chosen colormap from Matplotlib.
+"""
+function animate2d(files::Vector{String}, var; filedir="./", outdir="out/", vmin=-4, vmax=4,
+   cmap = matplotlib.cm.RdBu_r)
+   nfile = length(files)
+   x, y = let bd = load(filedir*files[1])
+      range(bd.x[1,1,1], bd.x[end,1,1], length=size(bd.x,1)),
+      range(bd.x[1,1,2], bd.x[1,end,2], length=size(bd.x,2))
+   end
+   # colormap norm
+   norm = matplotlib.colors.Normalize(vmin, vmax)
+   
+   fig = plt.figure(figsize=(8, 3), constrained_layout=true)
+   ax = plt.axes()
+
+   c = let
+      fakedata = zeros(Float32, length(x), length(y))
+      ax.pcolormesh(x, y, fakedata'; norm, cmap)
+   end
+   ax.set_xlabel("x [km]", fontsize=14)
+   ax.set_ylabel("y [km]", fontsize=14)
+   cb = colorbar(c; ax, orientation="horizontal", location="top", aspect=50)  
+   cb.ax.set_ylabel("Bz [nT]", fontsize=14)
+
+   for (i, file) in enumerate(files)
+      @info "$i in $nfile"
+      bd = load(joinpath(filedir, file))
+      c.set_array(bd[var]')
+   
+      title_str = @sprintf "t = %4.1f s" bd.head.time
+      ax.set_title(title_str)
+   
+      savefig(outdir*lpad(i, 4, '0')*".png", bbox_inches="tight", dpi=200)
+   end
+   
+   close()
+end
+
+"""
+    animate2d_with_streamline(files::Vector{String}, var::String, streamvars; kwargs...)
+
+Saving series of plots of `var`, together with streamplots of `streamvars`, from SWMF output
+`files`.
+
+# Keywords
+- `filedir::String="./"`: input SWMF file directory.
+- `outdir::String="out/"`: output directory.
+- `vmin=-4`: plot value lower bound.
+- `vmax=4`: plot value upper bound.
+- `cmap=matplotlib.cm.RdBu_r`: chosen colormap from Matplotlib.
+"""
+function animate2d_with_streamline(files::Vector{String}, var, streamvars="Bx;By";
+   filedir="./", outdir="out/", vmin=-4, vmax=4, cmap=matplotlib.cm.RdBu_r)
+   nfile = length(files)
+   x, y = let bd = load(filedir*files[1])
+      range(bd.x[1,1,1], bd.x[end,1,1], length=size(bd.x,1)),
+      range(bd.x[1,1,2], bd.x[1,end,2], length=size(bd.x,2))
+   end
+   # colormap norm
+   norm = matplotlib.colors.Normalize(vmin, vmax)
+
+   fig = plt.figure(figsize=(8, 3), constrained_layout=true)
+   ax = plt.axes()
+
+   @info "1 in $nfile"
+   c = ax.pcolormesh(x, y, bd[var]'; norm, cmap)
+
+   ax.set_xlabel("x [km]", fontsize=14)
+   ax.set_ylabel("y [km]", fontsize=14)
+   cb = colorbar(c; ax, orientation="horizontal", location="top", aspect=50)  
+   cb.ax.set_ylabel("Bz [nT]", fontsize=14)
+   title_str = @sprintf "t = %4.1f s" bd.head.time
+   ax.set_title(title_str)
+
+   st = streamplot(bd, streamvars, ax; color="gray", density=2)
+
+   save_and_clean!(1, outdir, st, ax)
+
+   for i in 2:nfile
+      @info "$i out of $nfile"
+      bd = load(joinpath(filedir, files[i]))
+      c.set_array(bd["Bz"]')
+
+      st = streamplot(bd, streamvars, ax; color="gray", density=2)
+      ax.set_xlim(x[1], x[end])
+      ax.set_ylim(y[1], y[end])
+
+      title_str = @sprintf "t = %4.1f s" bd.head.time
+      ax.set_title(title_str)
+   
+      save_and_clean!(i, outdir, st, ax)
+   end
+   
+   close()
+end
+
+function save_and_clean!(i::Int, outdir::String, st, ax)
+   savefig(outdir*lpad(i, 4, '0')*".png", bbox_inches="tight", dpi=200)
+   # Clean up streamlines
+   st.lines.remove()
+   for art in ax.get_children()
+      if !pybuiltin(:isinstance)(art, matplotlib.patches.FancyArrowPatch)
+         continue
+      end
+      art.remove()
+   end
+end
+
+####################
+# Data directory
+filedir = "./"
+# Plot variables
+var = "By"
+
+# Find simulation data
+files = let
+   pick_file = file -> startswith(file, "z") && endswith(file, ".out")
+   filter(pick_file, readdir(filedir))
+end
+
+animate2d(files, var)
+#animate2d_with_streamline(files, var, "Bx;By")
+```
diff --git a/docs/examples/config.json b/docs/examples/config.json
new file mode 100644
index 0000000..309f461
--- /dev/null
+++ b/docs/examples/config.json
@@ -0,0 +1,9 @@
+{
+   "theme": "nocoverlist",
+   "order": [
+     "basics"
+   ],
+   "properties": {
+     "notebook": "false"
+   }
+}
\ No newline at end of file
diff --git a/docs/examples/index.md b/docs/examples/index.md
new file mode 100644
index 0000000..8a38ca7
--- /dev/null
+++ b/docs/examples/index.md
@@ -0,0 +1,9 @@
+# [Examples](@id examples)
+
+This section contains thorough examples of using Batsrus.jl.
+
+---
+
+{{{democards}}}
+
+---
\ No newline at end of file
diff --git a/docs/make.jl b/docs/make.jl
index d3f72b7..ff8cd51 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -1,5 +1,13 @@
 using Batsrus, Documenter, PyPlot
 using Batsrus.UnitfulBatsrus
+using DemoCards
+
+branch = "master"
+# generate demo files
+demos, postprocess_cb, demo_assets = makedemos("examples"; branch)
+# if there are generated css assets, pass it to Documenter.HTML
+assets = String[]
+isnothing(demo_assets) || (push!(assets, demo_assets))
 
 makedocs(;
     modules=[Batsrus],
@@ -13,6 +21,7 @@ makedocs(;
     pages=[
         "Home" => "index.md",
         "Manual" => "man/manual.md",
+        "Examples" => demos,
         "Internal" => "man/internal.md",
         "Log" => "man/log.md",
     ],