diff --git a/examples/coremem.rs b/examples/coremem.rs
index 84d0b5e..b082e26 100644
--- a/examples/coremem.rs
+++ b/examples/coremem.rs
@@ -22,8 +22,6 @@ fn main() {
             // NB: different sources give pretty different values for this
             // NB: Simulation misbehaves for values > 10... Proably this model isn't so great.
             // Maybe use \eps or \xi instead of conductivity.
-            // No, I don't think that gives us what's needed. Rather, the calculation of static_ex
-            // is wrong (it's computed for the wrong value of time)
             state.get_mut(x, y).mat_mut().conductivity = (2.17).into();
         }
     }
diff --git a/src/lib.rs b/src/lib.rs
index 4a6a9a9..fddf3b8 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -287,11 +287,53 @@ pub trait Material {
 
 #[derive(Clone, Default)]
 pub struct GenericMaterial {
+    // Electrical parameters:
     pub conductivity: f64,
+    // Magnetic state:
+    /// Instantaneous magnetization in the z direction.
+    pub mz: f64,
+    // Ferromagnetic parameters:
+    /// Magnetic saturation (symmetric, so saturates to either +Ms or -Ms).
+    pub ms: f64,
+    /// Rate at which M changes w.r.t. H.
+    pub xi: f64,
 }
 
 impl Material for GenericMaterial {
     fn conductivity(&self) -> f64 {
         self.conductivity
     }
+    fn mz(&self) -> f64 {
+        self.mz
+    }
+    fn step_h(&mut self, context: &CellState, delta_bz: f64) -> f64 {
+        // COMSOL ferromagnetic model: https://www.comsol.com/blogs/accessing-external-material-models-for-magnetic-simulations/
+        //   delta(M) = xi * delta*(H)            # delta(M) = M(t+1) - M(t)
+        //   new_M = clamp(M+delta(M), -Ms, Ms))  # i.e. clamp to saturation
+        //   Where delta*(H) is H(t+1) - H(t), but clamping H(t) to <= 0 if positively saturated
+        //   and clamping H(t) to >= 0 if negatively saturated
+        let expected_delta_h = delta_bz / consts::MU0;
+        let expected_new_h = context.hz() + expected_delta_h;
+        // let unsaturated_new_m = self.mz + self.xi * unsaturated_delta_h;
+        let unsaturated_new_m = if self.mz == self.ms {
+            // positively saturated
+            let delta_h = expected_new_h.max(0.0) - context.hz().max(0.0);
+            self.mz + self.xi * delta_h
+        } else if self.mz == -self.ms {
+            // negatively saturated
+            let delta_h = expected_new_h.min(0.0) - context.hz().min(0.0);
+            self.mz + self.xi * delta_h
+        } else {
+            // not saturated
+            self.mz + self.xi * expected_delta_h
+        };
+        let new_m = unsaturated_new_m.min(self.ms).max(-self.ms);
+
+        let delta_m = new_m - self.mz;
+        self.mz = new_m;
+        // B = mu0*(H + M)
+        // \Delta B = mu0*(\Delta H + \Delta M)
+        let delta_h = delta_bz / consts::MU0 - delta_m;
+        context.hz() + delta_h
+    }
 }
diff --git a/src/render.rs b/src/render.rs
index 1a91ec8..f2b8b01 100644
--- a/src/render.rs
+++ b/src/render.rs
@@ -50,7 +50,8 @@ impl ColorTermRenderer {
             for x in 0..state.width() {
                 let cell = state.get(x, y);
                 //let r = norm_color(cell.bz() * consts::C);
-                let r = 0;
+                //let r = 0;
+                let r = norm_color(cell.mat.mz*50.0);
                 let b = (55.0*cell.mat().conductivity).min(255.0) as u8;
                 //let b = 0;
                 //let g = norm_color(cell.ex());