diff --git a/meanas/fdmath/__init__.py b/meanas/fdmath/__init__.py
index 76bd908..cb47c79 100644
--- a/meanas/fdmath/__init__.py
+++ b/meanas/fdmath/__init__.py
@@ -288,9 +288,9 @@ If we discretize both space (m,n,p) and time (l), Maxwell's equations become
 
  $$ \\begin{align*}
   \\tilde{\\nabla} \\times \\tilde{E}_{l,\\vec{r}} &= -\\tilde{\\partial}_t \\hat{B}_{l-\\frac{1}{2}, \\vec{r} + \\frac{1}{2}}
-                                                                         + \\hat{M}_{l-1, \\vec{r} + \\frac{1}{2}}  \\\\
-  \\hat{\\nabla} \\times \\hat{H}_{l,\\vec{r} + \\frac{1}{2}} &= \\hat{\\partial}_t \\tilde{D}_{l, \\vec{r}}
-                                                                   + \\tilde{J}_{l-\\frac{1}{2},\\vec{r}} \\\\
+                                                                         - \\hat{M}_{l, \\vec{r} + \\frac{1}{2}}  \\\\
+  \\hat{\\nabla} \\times \\hat{H}_{l-\\frac{1}{2},\\vec{r} + \\frac{1}{2}} &= \\hat{\\partial}_t \\tilde{D}_{l, \\vec{r}}
+                                                                             + \\tilde{J}_{l-\\frac{1}{2},\\vec{r}} \\\\
   \\tilde{\\nabla} \\cdot \\hat{B}_{l-\\frac{1}{2}, \\vec{r} + \\frac{1}{2}} &= 0 \\\\
   \\hat{\\nabla} \\cdot \\tilde{D}_{l,\\vec{r}} &= \\rho_{l,\\vec{r}}
  \\end{align*} $$
@@ -311,9 +311,9 @@ and \\( \\epsilon \\) and \\( \\mu \\) are the dielectric permittivity and magne
 The above is Yee's algorithm, written in a form analogous to Maxwell's equations.
 The time derivatives can be expanded to form the update equations:
 
-    [code: Maxwell's equations]
-    H[i, j, k] -= (curl_forward(E[t])[i, j, k] - M[t, i, j, k]) /      mu[i, j, k]
-    E[i, j, k] += (curl_back(   H[t])[i, j, k] + J[t, i, j, k]) / epsilon[i, j, k]
+    [code: Maxwell's equations updates]
+    H[i, j, k] -= dt * (curl_forward(E)[i, j, k] + M[t, i, j, k]) /      mu[i, j, k]
+    E[i, j, k] += dt * (curl_back(   H)[i, j, k] + J[t, i, j, k]) / epsilon[i, j, k]
 
 Note that the E-field fore-vector and H-field back-vector are offset by a half-cell, resulting
 in distinct locations for all six E- and H-field components:
@@ -383,7 +383,7 @@ $$
   \\begin{align*}
   \\tilde{\\nabla} \\times \\tilde{E}_{l,\\vec{r}} &=
       -\\tilde{\\partial}_t \\hat{B}_{l-\\frac{1}{2}, \\vec{r} + \\frac{1}{2}}
-                          + \\hat{M}_{l-1, \\vec{r} + \\frac{1}{2}}  \\\\
+                          - \\hat{M}_{l-1, \\vec{r} + \\frac{1}{2}}  \\\\
   \\mu^{-1}_{\\vec{r} + \\frac{1}{2}} \\cdot \\tilde{\\nabla} \\times \\tilde{E}_{l,\\vec{r}} &=
     -\\tilde{\\partial}_t \\hat{H}_{l-\\frac{1}{2}, \\vec{r} + \\frac{1}{2}}  \\\\
   \\hat{\\nabla} \\times (\\mu^{-1}_{\\vec{r} + \\frac{1}{2}} \\cdot \\tilde{\\nabla} \\times \\tilde{E}_{l,\\vec{r}}) &=
@@ -488,11 +488,16 @@ $$
   K_x^2 + K_y^2 + K_z^2 = \\Omega^2 \\mu \\epsilon = \\Omega^2 / c^2
 $$
 
+where \\( c = \\sqrt{\\mu \\epsilon} \\).
+
 Assuming real \\( (k_x, k_y, k_z), \\omega \\) will be real only if
 
 $$ c^2 \\Delta_t^2 = \\frac{\\Delta_t^2}{\\mu \\epsilon} < 1/(\\frac{1}{\\Delta_x^2} + \\frac{1}{\\Delta_y^2} + \\frac{1}{\\Delta_z^2}) $$
 
-If \\( \\Delta_x = \\Delta_y = \\Delta_z \\), this simplifies to \\( c \\Delta_t  < \\frac{\\Delta_x}{\\sqrt{3}} \\).
+If \\( \\Delta_x = \\Delta_y = \\Delta_z \\), this simplifies to \\( c \\Delta_t < \\Delta_x / \\sqrt{3} \\).
+This last form can be interpreted as enforcing causality; the distance that light
+travels in one timestep (i.e., \\( c \\Delta_t \\)) must be less than the diagonal
+of the smallest cell ( \\( \\Delta_x / \\sqrt{3} \\) when on a uniform cubic grid).
 
 
 Grid description
diff --git a/meanas/fdtd/__init__.py b/meanas/fdtd/__init__.py
index 2d1588e..3d498a9 100644
--- a/meanas/fdtd/__init__.py
+++ b/meanas/fdtd/__init__.py
@@ -1,6 +1,9 @@
 """
 Utilities for running finite-difference time-domain (FDTD) simulations
 
+See the discussion of `Maxwell's Equations` in `meanas.fdmath` for basic
+mathematical background.
+
 
 Timestep
 ========
@@ -19,13 +22,120 @@ Based on this, we can set
 The `dx_min`, `dy_min`, `dz_min` should be the minimum value across both the base and derived grids.
 
 
-Poynting Vector
-===============
-# TODO
+Poynting Vector and Energy Conservation
+=======================================
+
+Let
+
+$$ \\begin{align*}
+  \\tilde{S}_{l, l', \\vec{r}} &=& &\\tilde{E}_{l, \\vec{r}} \\otimes \\hat{H}_{l', \\vec{r} + \\frac{1}{2}}  \\\\
+  &=&  &\\vec{x} (\\tilde{E}^y_{l,m+1,n,p} \\hat{H}^z_{l',\\vec{r} + \\frac{1}{2}} - \\tilde{E}^z_{l,m+1,n,p} \\hat{H}^y_{l', \\vec{r} + \\frac{1}{2}}) \\\\
+  & &+ &\\vec{y} (\\tilde{E}^z_{l,m,n+1,p} \\hat{H}^x_{l',\\vec{r} + \\frac{1}{2}} - \\tilde{E}^x_{l,m,n+1,p} \\hat{H}^z_{l', \\vec{r} + \\frac{1}{2}}) \\\\
+  & &+ &\\vec{z} (\\tilde{E}^x_{l,m,n,p+1} \\hat{H}^y_{l',\\vec{r} + \\frac{1}{2}} - \\tilde{E}^y_{l,m,n,p+1} \\hat{H}^z_{l', \\vec{r} + \\frac{1}{2}})
+   \\end{align*}
+$$
+
+where \\( \\vec{r} = (m, n, p) \\) and \\( \\otimes \\) is a modified cross product
+in which the \\( \\tilde{E} \\) terms are shifted as indicated.
+
+By taking the divergence and rearranging terms, we can show that
+
+$$
+  \\begin{align*}
+  \\hat{\\nabla} \\cdot \\tilde{S}_{l, l', \\vec{r}}
+   &= \\hat{\\nabla} \\cdot (\\tilde{E}_{l, \\vec{r}} \\otimes \\hat{H}_{l', \\vec{r} + \\frac{1}{2}})  \\\\
+   &= \\hat{H}_{l', \\vec{r} + \\frac{1}{2}} \\cdot \\tilde{\\nabla} \\times \\tilde{E}_{l, \\vec{r}} -
+      \\tilde{E}_{l, \\vec{r}} \\cdot \\hat{\\nabla} \\times \\hat{H}_{l', \\vec{r} + \\frac{1}{2}} \\\\
+   &= \\hat{H}_{l', \\vec{r} + \\frac{1}{2}} \\cdot
+            (-\\tilde{\\partial}_t \\mu_{\\vec{r} + \\frac{1}{2}} \\hat{H}_{l - \\frac{1}{2}, \\vec{r} + \\frac{1}{2}} -
+                \\hat{M}_{l-1, \\vec{r} + \\frac{1}{2}}) -
+      \\tilde{E}_{l, \\vec{r}} \\cdot (\\hat{\\partial}_t \\tilde{\\epsilon}_\\vec{r} \\tilde{E}_{l'+\\frac{1}{2}, \\vec{r}} +
+                \\tilde{J}_{l', \\vec{r}}) \\\\
+   &= \\hat{H}_{l'} \\cdot (-\\mu / \\Delta_t)(\\hat{H}_{l + \\frac{1}{2}} - \\hat{H}_{l - \\frac{1}{2}}) -
+      \\tilde{E}_l \\cdot (\\epsilon / \\Delta_t )(\\tilde{E}_{l'+\\frac{1}{2}} - \\tilde{E}_{l'-\\frac{1}{2}})
+      - \\hat{H}_{l'} \\cdot \\hat{M}_{l-1} - \\tilde{E}_l \\cdot \\tilde{J}_{l'} \\\\
+  \\end{align*}
+$$
+
+where in the last line the spatial subscripts have been dropped to emphasize
+the time subscripts \\( l, l' \\), i.e.
+
+$$
+  \\begin{align*}
+  \\tilde{E}_l &= \\tilde{E}_{l, \\vec{r}} \\\\
+  \\hat{H}_l &= \\tilde{H}_{l, \\vec{r} + \\frac{1}{2}}  \\\\
+  \\tilde{\\epsilon} &= \\tilde{\\epsilon}_\\vec{r}  \\\\
+  \\end{align*}
+$$
+
+etc.
+For \\( l' = l + \\frac{1}{2} \\) we get
+
+$$
+  \\begin{align*}
+  \\hat{\\nabla} \\cdot \\tilde{S}_{l, l + \\frac{1}{2}}
+   &= \\hat{H}_{l + \\frac{1}{2}} \\cdot
+            (-\\mu / \\Delta_t)(\\hat{H}_{l + \\frac{1}{2}} - \\hat{H}_{l - \\frac{1}{2}}) -
+      \\tilde{E}_l \\cdot (\\epsilon / \\Delta_t)(\\tilde{E}_{l+1} - \\tilde{E}_l)
+      - \\hat{H}_{l'} \\cdot \\hat{M}_l - \\tilde{E}_l \\cdot \\tilde{J}_{l + \\frac{1}{2}} \\\\
+   &= (-\\mu / \\Delta_t)(\\hat{H}^2_{l + \\frac{1}{2}} - \\hat{H}_{l + \\frac{1}{2}} \\cdot \\hat{H}_{l - \\frac{1}{2}}) -
+      (\\epsilon / \\Delta_t)(\\tilde{E}_{l+1} \\cdot \\tilde{E}_l - \\tilde{E}^2_l)
+      - \\hat{H}_{l'} \\cdot \\hat{M}_l - \\tilde{E}_l \\cdot \\tilde{J}_{l + \\frac{1}{2}} \\\\
+   &= -(\\mu \\hat{H}^2_{l + \\frac{1}{2}}
+       +\\epsilon \\tilde{E}_{l+1} \\cdot \\tilde{E}_l) / \\Delta_t \\ \\
+      +(\\mu \\hat{H}_{l + \\frac{1}{2}} \\cdot \\hat{H}_{l - \\frac{1}{2}}
+       +\\epsilon \\tilde{E}^2_l) / \\Delta_t \\ \\
+      - \\hat{H}_{l+\\frac{1}{2}} \\cdot \\hat{M}_l \\ \\
+      - \\tilde{E}_l \\cdot \\tilde{J}_{l+\\frac{1}{2}} \\\\
+  \\end{align*}
+$$
+
+and for \\( l' = l - \\frac{1}{2} \\),
+
+$$
+  \\begin{align*}
+  \\hat{\\nabla} \\cdot \\tilde{S}_{l, l - \\frac{1}{2}}
+   &=  (\\mu \\hat{H}^2_{l - \\frac{1}{2}}
+       +\\epsilon \\tilde{E}_{l-1} \\cdot \\tilde{E}_l) / \\Delta_t \\ \\
+      -(\\mu \\hat{H}_{l + \\frac{1}{2}} \\cdot \\hat{H}_{l - \\frac{1}{2}}
+       +\\epsilon \\tilde{E}^2_l) / \\Delta_t \\ \\
+      - \\hat{H}_{l-\\frac{1}{2}} \\cdot \\hat{M}_l \\ \\
+      - \\tilde{E}_l \\cdot \\tilde{J}_{l-\\frac{1}{2}} \\\\
+  \\end{align*}
+$$
+
+These two results form the discrete time-domain analogue to Poynting's theorem.
+They hint at the expressions for the energy, which can be calculated at the same
+time-index as either the E or H field:
+
+$$
+ \\begin{align*}
+ U_l &= \\epsilon \\tilde{E}^2_l + \\mu \\hat{H}_{l + \\frac{1}{2}} \\cdot \\hat{H}_{l - \\frac{1}{2}} \\\\
+ U_{l + \\frac{1}{2}} &= \\epsilon \\tilde{E}_l \\cdot \\tilde{E}_{l + 1} + \\mu \\hat{H}^2_{l + \\frac{1}{2}} \\\\
+ \\end{align*}
+$$
+
+Rewriting the Poynting theorem in terms of the energy expressions,
+
+$$
+  \\begin{align*}
+  (U_{l+\\frac{1}{2}} - U_l) / \\Delta_t
+   &= -\\hat{\\nabla} \\cdot \\tilde{S}_{l, l + \\frac{1}{2}} \\ \\
+      - \\hat{H}_{l+\\frac{1}{2}} \\cdot \\hat{M}_l \\ \\
+      - \\tilde{E}_l \\cdot \\tilde{J}_{l+\\frac{1}{2}} \\\\
+  (U_l - U_{l-\\frac{1}{2}}) / \\Delta_t
+   &= -\\hat{\\nabla} \\cdot \\tilde{S}_{l, l - \\frac{1}{2}} \\ \\
+      - \\hat{H}_{l-\\frac{1}{2}} \\cdot \\hat{M}_l \\ \\
+      - \\tilde{E}_l \\cdot \\tilde{J}_{l-\\frac{1}{2}} \\\\
+ \\end{align*}
+$$
+
+This result is exact an should practically hold to within numerical precision. No time-
+or spatial-averaging is necessary.
+
+Note that each value of \\( J \\) contributes to the energy twice (i.e. once per field update)
+despite only causing the value of \\( E \\) to change once (same for \\( M \\) and \\( H \\)).
 
-Energy conservation
-===================
-# TODO
 
 Boundary conditions
 ===================
diff --git a/meanas/fdtd/base.py b/meanas/fdtd/base.py
index fd21642..e16a53f 100644
--- a/meanas/fdtd/base.py
+++ b/meanas/fdtd/base.py
@@ -73,9 +73,9 @@ def maxwell_h(dt: float, dxes: dx_lists_t = None) -> fdfield_updater_t:
 
     The full update should be
 
-        H -= (curl_forward(E[t]) - M) / mu
+        H -= (curl_forward(E[t]) + M) / mu
 
-    which requires an additional step of `H += M / mu` which is not performed
+    which requires an additional step of `H -= M / mu` which is not performed
     by the generated function; this step can be omitted if there is no magnetic
     current `M`.