diff --git a/Project.toml b/Project.toml
index 2016bdc2..1cd7edfb 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,21 +5,26 @@ version = "0.0.0-DEV"
 
 [deps]
 Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"
+Bonito = "824d6782-a2ef-11e9-3a09-e5662e0c26f8"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
 Gridap = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
 GridapGmsh = "3025c34a-b394-11e9-2a55-3fee550c04c8"
 GridapMakie = "41f30b06-6382-4b60-a5f7-79d86b35bf5d"
 GridapPETSc = "bcdc36c2-0c3e-11ea-095a-c9dadae499f1"
+Gumbo = "708ec375-b3d6-5a57-a7ce-8257bf98657a"
+HTTP = "cd3eb016-35fb-5094-929b-558a96fad6f3"
 IJulia = "7073ff75-c697-5162-941a-fcdaad2a7d2a"
 JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899"
 LsqFit = "2fda8390-95c7-5789-9bda-21331edee243"
 MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
 MUMPS = "55d2b088-9f4e-11e9-26c0-150b02ea6a46"
+Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 PhysicalConstants = "5ad8b20f-a522-5ce9-bfc9-ddf1d5bda6ab"
 Pickle = "fbb45041-c46e-462f-888f-7c521cafbc2c"
 PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
+WGLMakie = "276b4fcb-3e11-5398-bf8b-a0c2d153d008"
 
 [compat]
 Arpack = "=0.5.3"
diff --git a/docs/julia/heater_3d.jl b/docs/julia/heater_3d.jl
index b51087ba..dcdb8d0b 100644
--- a/docs/julia/heater_3d.jl
+++ b/docs/julia/heater_3d.jl
@@ -18,9 +18,12 @@
 # # 3D thermal phase shifter
 
 # %% tags=["remove-stderr", "hide-input", "hide-output"]
-using CairoMakie
+using GridapMakie
 using Printf
 
+include("./js.jl")
+using .JS: FigureContainer
+
 using Gridap
 using GridapGmsh
 using Gridap.Geometry
@@ -29,8 +32,201 @@ using GridapPETSc
 using Femwell.Maxwell.Electrostatic
 using Femwell.Thermal
 
+using IJulia
+using CairoMakie
+using Markdown
+
+using WGLMakie
+WGLMakie.activate!()
+
+Makie.inline!(true) # Make sure to inline plots into Documenter output!
+
 dir = @__DIR__
-run(`python $dir/heater_3d_mesh.py`)
+# run(`python $dir/heater_3d_mesh.py`)
+
+# %% tags=["remove-stderr", "hide-input", "hide-output"]
+
+function visualize_mesh(model)
+    volume_tags = [
+        "box",
+        "core",
+        "heater",
+        "via2",
+        "via1",
+        "metal3",
+        "metal2",
+        "metal3#e1",
+        "metal3#e2",
+    ]
+    colors = [:grey, :black, :orange, :green, :green, :brown, :brown, :blue, :blue]
+    alphas = [1.0, 1.0, 0.2, 0.5, 0.5, 0.1, 0.1, 0.75, 0.75]
+    fig = Figure(fontsize = 16)
+    xlims = (-20e-6, 60e-6)
+    ylims = (-5e-6, 5e-6)
+    zlims = (-5e-6, 5e-6)
+    aspect = (xlims[2] - xlims[1]) / (ylims[2] - ylims[1])
+    tickfunc = values -> (x -> string(x * 1e6)).(values)
+
+    # Mesh
+    ax =
+        fig[1, 2] = Axis3(
+            fig[1, 2],
+            aspect = (aspect, 1.0, 1.0),
+            xlabel = "x (μm)",
+            ylabel = "\ny (μm)",
+            zlabel = "z (μm)\n\n",
+            title = "Heater Mesh",
+            limits = (xlims, ylims, zlims),
+            titlesize = 28,
+            xtickformat = tickfunc,
+            ytickformat = tickfunc,
+            ztickformat = tickfunc,
+            azimuth = 51 / 40 * pi + pi / 20,
+            elevation = pi / 8 + pi / 10,
+        )
+    Camera3D(ax.scene)
+
+    for (i, tag) in enumerate(volume_tags)
+        ∂Ω_region = BoundaryTriangulation(model, tags = [tag])
+        wireframe!(
+            fig[1, 2],
+            ∂Ω_region,
+            color = colors[i],
+            alpha = alphas[i],
+            linewidth = 0.1,
+            fontsize = 16,
+        )
+    end
+    wf = [PolyElement(color = color) for color in colors]
+    Legend(
+        fig[1, 1],
+        wf,
+        volume_tags,
+        halign = :left,
+        valign = :top,
+        framevisible = false,
+        labelsize = 16,
+    )
+
+    display(FigureContainer(fig))
+end
+
+function visualize_field(model, field, label; plotargs...)
+    volume_tags = [
+        "box",
+        "core",
+        "heater",
+        "via2",
+        "via1",
+        "metal3",
+        "metal2",
+        "metal3#e1",
+        "metal3#e2",
+    ]
+    fontsize = 10
+    fig = Figure(fontsize = fontsize, resolution = (800, 500))
+    xlims = (-20e-6, 60e-6)
+    ylims = (-5e-6, 5e-6)
+    zlims = (-5e-6, 5e-6)
+    aspect = (xlims[2] - xlims[1]) / (ylims[2] - ylims[1])
+    tickfunc = values -> (x -> string(x * 1e6)).(values)
+
+    # Setup plot
+    Ω = Triangulation(model, tags = volume_tags)
+    ax =
+        fig[1, 1] = Axis3(
+            fig[1, 1],
+            aspect = (aspect, 1.0, 1.0),
+            xlabel = "x (μm)",
+            ylabel = "y (μm)",
+            zlabel = "z (μm)",
+            title = label,
+            limits = (xlims, ylims, zlims),
+            titlesize = 18,
+            xtickformat = tickfunc,
+            ytickformat = tickfunc,
+            ztickformat = tickfunc,
+            azimuth = 51 / 40 * pi + pi / 20,
+            elevation = pi / 8 + pi / 10,
+        )
+    Camera3D(ax.scene)
+
+    # Plot
+    plt = plot!(fig[1, 1], Ω, field, shading = true; plotargs...)
+    Colorbar(fig[1, 2], plt, label = label, vertical = true, flipaxis = false)
+
+    display(FigureContainer(fig))
+end
+
+function plot_time(uₕₜ::Gridap.ODEs.TransientFETools.TransientFESolution; label = "")
+    # Set up figure
+    volume_tags = [
+        "box",
+        "core",
+        "heater",
+        "via2",
+        "via1",
+        "metal3",
+        "metal2",
+        "metal3#e1",
+        "metal3#e2",
+    ]
+    fig = CairoMakie.Figure()
+    xlims = (-20e-6, 60e-6)
+    ylims = (-5e-6, 5e-6)
+    zlims = (-5e-6, 5e-6)
+    aspect = (xlims[2] - xlims[1]) / (ylims[2] - ylims[1])
+    tickfunc = values -> (x -> string(x * 1e6)).(values)
+    Ω = Triangulation(model, tags = volume_tags)
+
+    # Set up axis
+    it = iterate(uₕₜ)
+    uh = Any[nothing]
+    time = Any[nothing]
+    state = Any[nothing]
+    (uh[1], time[1]), state[1] = it
+    t_uh_time = Observable(time[1])
+    timetitle = lift(time -> label * @sprintf("   t = %.2f ms", time * 1e3), t_uh_time)
+    Axis3(
+        fig[1, 1],
+        aspect = (aspect, 1.0, 1.0),
+        xlabel = "x (μm)",
+        ylabel = "\ny (μm)",
+        zlabel = "z (μm)",
+        title = timetitle,
+        limits = (xlims, ylims, zlims),
+        xtickformat = tickfunc,
+        ytickformat = tickfunc,
+        ztickformat = tickfunc,
+        azimuth = 51 / 40 * pi + pi / 20,
+        elevation = pi / 8 + pi / 10,
+    )
+
+    # Plot new iteration
+    u = lift(t_uh_time) do t
+        while (time[1] < t) || (time[1] ≈ t)
+            if time[1] ≈ t
+                return uh[1]
+            end
+            it = iterate(uₕₜ, state[1])
+            if isnothing(it)
+                return uh[1]
+            end
+            (uh[1], time[1]), state[1] = it
+        end
+        return uh[1]
+    end
+
+    # Create the actual plot
+    plt = plot!(fig[1, 1], Ω, u, colorrange = (0, 110.0), colormap = Reverse(:RdBu))
+    Colorbar(fig[1, 2], plt, vertical = true, label = "Temperature (K)")
+    colsize!(fig.layout, 1, Aspect(1, 1.0))
+    t_uh_time, fig
+end
+
+# %% tags=["remove-stderr"]
+model = GmshDiscreteModel("$dir/mesh.msh")
+visualize_mesh(model)
 
 # %% [markdown]
 # Let's start with defining the constants for the simulation:
@@ -77,7 +273,6 @@ thermal_diffisitivities = [
 # and create functions which work on the tags
 
 # %% tags=["remove-stderr", "hide-output"]
-model = GmshDiscreteModel("mesh.msh")
 Ω = Triangulation(model)
 dΩ = Measure(Ω, 1)
 
@@ -134,13 +329,53 @@ println(
     " K",
 )
 
+# %% tags=["remove-stderr"]
+visualize_field(
+    model,
+    potential(p0),
+    "Potential";
+    colormap = :viridis,
+    colorrange = (0, 0.5),
+    fontsize = 14,
+)
+
+# %% tags=["remove-stderr"]
+visualize_field(
+    model,
+    current_density(p0),
+    "Current Density";
+    colormap = :viridis,
+    colorrange = (0, 1.75e10),
+    fontsize = 14,
+)
+
+# %% tags=["remove-stderr"]
+visualize_field(
+    model,
+    temperature(T0),
+    "Temperature";
+    colormap = Reverse(:RdBu),
+    colorrange = (0, 120.0),
+    fontsize = 14,
+)
+
+# %% tags=["remove-stderr"]
+visualize_field(
+    model,
+    heat_flux(T0),
+    "Heat Flux";
+    colormap = Reverse(:RdBu),
+    colorrange = (3.0e4, 4e9),
+    fontsize = 14,
+)
+
 # %% [markdown]
-# And we write the fields to a file for visualisation using paraview:
+# Alternatively, we write the fields to a file for visualisation using paraview:
 
 # %% tags=["remove-stderr", "hide-output"]
 writevtk(
     Ω,
-    "results",
+    "$dir/results",
     cellfields = [
         "potential" => potential(p0),
         "current" => current_density(p0),
@@ -157,6 +392,7 @@ writevtk(
 # For the cooldown simulation we start with the steady state temperature profile and no heating.
 
 # %% tags=["remove-stderr"]
+CairoMakie.activate!()
 figure = Figure()
 ax = Axis(
     figure[1, 1],
@@ -164,6 +400,17 @@ ax = Axis(
     xlabel = "Time / ms",
 )
 
+# %% [markdown]
+# Here we set the parameters for the transient simulation
+# %% tags=["remove-stderr", "hide-output"]
+Δt = 2e-7
+t_end = 2e-4
+plot_every = 200
+timestamps = collect(Δt:Δt*plot_every:t_end)
+total_render_time = 5.0
+fps = ceil(Int, length(timestamps) / total_render_time)
+html_sources = String[]
+lins = CairoMakie.Lines[]
 for (label, power_factor, temperature_factor) in [("heatup", 1, 0), ("cooldown", 0, 1)]
     uₕₜ = calculate_temperature_transient(
         ϵ_conductivities ∘ τ,
@@ -171,10 +418,17 @@ for (label, power_factor, temperature_factor) in [("heatup", 1, 0), ("cooldown",
         power_density(p0) * power_factor,
         boundary_temperatures,
         temperature(T0) * temperature_factor,
-        2e-7,
-        2e-4,
+        Δt,
+        t_end,
     )
 
+    # Create animation
+    t_uh_time, fig = plot_time(uₕₜ; label)
+    record(fig, "$dir/animation-$label.gif", timestamps, framerate = fps) do this_t
+        t_uh_time[] = this_t
+    end
+
+    # Optionally write the solution to a file
     #createpvd("poisson_transient_solution_$label") do pvd
     #    for (uₕ, t) in uₕₜ
     #        pvd[t] = createvtk(
@@ -187,10 +441,25 @@ for (label, power_factor, temperature_factor) in [("heatup", 1, 0), ("cooldown",
 
     sums = [(t, ∑(∫(u)dΩ_w) / silicon_volume) for (u, t) in uₕₜ]
     t, s = getindex.(sums, 1), getindex.(sums, 2)
-    lines!(ax, t * 1e3, s, label = label)
+    push!(lins, lines!(ax, t * 1e3, s, label = label))
 end
 
-axislegend(position = :rc)
+# %% [markdown]
+# If we display the animations we can visualize the heatup and cooldown of the device:
+# ![heatup animation](animation-heatup.gif)
+# ![cooldown animation](animation-cooldown.gif)
+
+# %% [markdown]
+# Finally, we can plot an average of the temperature over time: 
+# %% tags=["remove-stderr"]
+Legend(
+    figure[1, 1],
+    lins,
+    ["heatup", "cooldown"],
+    halign = :right,
+    valign = :center,
+    labelsize = 16,
+)
 display(figure)
 
 
@@ -203,7 +472,7 @@ display(figure)
 #   end
 # ```
 # 
-# And then add to each of the calculate-calls the sovler parameter:
+# And then add to each of the calculate-calls the solver parameter:
 #
 # ```julia
 #   solver = PETScLinearSolver()
diff --git a/docs/julia/js.jl b/docs/julia/js.jl
new file mode 100644
index 00000000..7b0bf675
--- /dev/null
+++ b/docs/julia/js.jl
@@ -0,0 +1,77 @@
+module JS
+
+using Gumbo, HTTP, Bonito
+using Makie, WGLMakie
+export FigureContainer
+
+abstract type BackendType end
+struct WGLMakieType <: BackendType end
+
+function backend_type(::Module)
+    if Makie.current_backend() === WGLMakie
+        return WGLMakieType
+    else
+        return BackendType
+    end
+end
+
+struct FigureContainer{T<:BackendType}
+    fig::Figure
+end
+FigureContainer(fig::Figure) = FigureContainer{backend_type(Makie.current_backend())}(fig)
+
+js_counter = 1
+dir = @__DIR__
+
+function separate_javascript!(element::HTMLElement)
+    for (i, child) in enumerate(element.children)
+        if (typeof(child) <: HTMLElement{:script}) &&
+           haskey(child.attributes, "src") &&
+           contains(child.attributes["src"], "data:application/javascript;base64,")
+            src_raw = replace(
+                child.attributes["src"],
+                "data:application/javascript;base64," => "",
+            )
+            js_content = HTTP.base64decode(src_raw)
+            global js_counter += 1
+            js_filename = "./js/js-$(js_counter).js"
+            js_write = "$dir/../_build/html/julia/js/js-$(js_counter).js"
+            write(js_write, js_content)
+            element.children[i] = HTMLElement{:script}(
+                Vector{HTMLNode}[],
+                element,
+                Dict{String,String}("src" => js_filename, "type" => "module"),
+            )
+        end
+        if typeof(child) <: HTMLElement
+            separate_javascript!(child)
+        end
+    end
+    return element
+end
+
+function separate_javascript(html_source::AbstractString)
+    if !isdir("$dir/../_build/html/")
+        mkdir("$dir/../_build/html/")
+    end
+    if !isdir("$dir/../_build/html/julia/")
+        mkdir("$dir/../_build/html/julia/")
+    end
+    if !isdir("$dir/../_build/html/julia/js/")
+        mkdir("$dir/../_build/html/julia/js/")
+    end
+    parsed_html = parsehtml(html_source)
+    div_element = separate_javascript!(parsed_html.root.children[2].children[1])
+end
+
+function Base.display(fig::FigureContainer{<:BackendType})
+    display(fig.fig)
+end
+
+function Base.display(fig::FigureContainer{WGLMakieType})
+    html_str = string(DOM.div(fig.fig))
+    html_source_new = separate_javascript(html_str)
+    display("text/html", string(html_source_new))
+end
+
+end