diff --git a/docs/julia/heater_3d.jl b/docs/julia/heater_3d.jl
index 4f8de690..f46f7ca4 100644
--- a/docs/julia/heater_3d.jl
+++ b/docs/julia/heater_3d.jl
@@ -29,7 +29,7 @@ using Femwell.Maxwell.Electrostatic
 using Femwell.Thermal
 
 dir = @__DIR__
-read(`python $dir/heater_3d_mesh.py`)
+#read(`python $dir/heater_3d_mesh.py`)
 
 
 model = GmshDiscreteModel("mesh.msh")
@@ -143,7 +143,11 @@ GridapPETSc.with(args = split(options)) do
 
     Ω_w = Triangulation(model, tags = "core")
     dΩ_w = Measure(Ω_w, 1)
-    sums = [(t, ∑(∫(u)dΩ_w)) for (u, t) in uₕₜ]
-    lines(sums)
+    sums = [(t, ∑(∫(u)dΩ_w) / ∑(∫(1)dΩ_w)) for (u, t) in uₕₜ]
 
+    figure = Figure()
+    ax = Axis(figure[1, 1], ylabel = "Temperature / K", xlabel = "time / ms")
+
+    t, s = getindex.(sums, 1), getindex.(sums, 2)
+    lines!(ax, t * 1e3, s)
 end
diff --git a/src/Thermal/Thermal.jl b/src/Thermal/Thermal.jl
index cc1ea70f..a598d6a0 100644
--- a/src/Thermal/Thermal.jl
+++ b/src/Thermal/Thermal.jl
@@ -89,8 +89,8 @@ function calculate_temperature_transient(
         b₀(t, v) = ∫(power_density(t) * v)dΩ
         op_C = TransientConstantMatrixFEOperator(m₀, a₀, b₀, U, V)
     else
-        b₀(v) = ∫(power_density * v)dΩ
-        op_C = TransientConstantFEOperator(m₀, a₀, b₀, U, V)
+        b(v) = ∫(power_density * v)dΩ
+        op_C = TransientConstantFEOperator(m₀, a₀, b, U, V)
     end
 
     θ = 0.5