diff --git a/physica-manual.pdf b/physica-manual.pdf
index 90ae372..6f0ba24 100644
Binary files a/physica-manual.pdf and b/physica-manual.pdf differ
diff --git a/physica-manual.typ b/physica-manual.typ
index 0be2ee7..efeb21b 100644
--- a/physica-manual.typ
+++ b/physica-manual.typ
@@ -852,6 +852,51 @@ $ isotope("Bi",a:211,z:83) --> isotope("Tl",a:207,z:81) + isotope("He",a:4,z:2)
 *(4)* #hl(`isotope("Tl",a:207,z:81) --> isotope("Pb",a:207,z:82) + isotope(e,a:0,z:-1)`)
 $ isotope("Tl",a:207,z:81) --> isotope("Pb",a:207,z:82) + isotope(e,a:0,z:-1) $
 
+=== The n-th term in Taylor series
+
+#v(1em)
+
+Function: `taylorterm(`_func_, _x_, _x0_, _idx_`)`.
+- _func_: the function e.g. `f`, `(f+g)`,
+- _x_: the variable name e.g. `x`,
+- _x0_: the variable value at the expansion point e.g. `x_0`, `(1+a)`,
+- _idx_: the index of the term, e.g. `0`, `1`, `2`, `n`, `(n+1)`.
+If _x0_ or _idx_ is in a parenthesis e.g. `(1+k)`, then the function
+automatically removes the outer parenthesis where appropriate.
+
+#align(center, [*Examples*])
+
+#grid(
+  columns: (50%, 50%),
+  row-gutter: 1em,
+  column-gutter: 2em,
+  [
+    *(1)* #hl(`taylorterm(f,x,x_0,0)`) \
+    $ taylorterm(f,x,x_0,0) $
+  ],
+  [
+    *(2)* #hl(`taylorterm(f,x,x_0,1)`) \
+    $ taylorterm(f,x,x_0,1) $
+  ],
+  [
+    *(3)* #hl(`taylorterm(f,x,(1+a),2)`) \
+    $ taylorterm(f,x,(1+a),2) $
+  ],
+  [
+    *(4)* #hl(`taylorterm(f,x,x_0,n)`) \
+    $ taylorterm(f,x,x_0,n) $
+  ],
+  [
+    *(5)* #hl(`taylorterm(F,x^nu,x^nu_0,n)`) \
+    $ taylorterm(F,x^nu,x^nu_0,n) $
+  ],
+  [
+    *(6)* #hl(`taylorterm(f_p,x,x_0,(n+1))`) \
+    $ taylorterm(f_p,x,x_0,(n+1)) $
+  ],
+)
+
+
 === Signal sequences (digital timing diagrams)
 
 In engineering, people often need to draw digital timing diagrams for signals, like $signals("1|0|1|0")$.
diff --git a/physica.typ b/physica.typ
index b4a8271..9b742a6 100644
--- a/physica.typ
+++ b/physica.typ
@@ -288,8 +288,8 @@
 #let jmat = jacobianmatrix
 
 #let hessianmatrix(fs, xs, delim:"(") = {
-  assert(type(fs) == array, message: "expecting a one-element array")
-  assert(fs.len() == 1, message: "expecting only one function name")
+  assert(type(fs) == array, message: "usage: hessianmatrix(f; x, y...)")
+  assert(fs.len() == 1, message: "usage: hessianmatrix(f; x, y...)")
   let f = fs.at(0)
   assert(type(xs) == array, message: "expecting an array of variable names")
   let row_arrays = ()  // array of arrays
@@ -663,6 +663,25 @@
   math.attach((T,hphantom(sym.zwj)).join(), t: uppers.join(), b: lowers.join())
 }
 
+#let taylorterm(fn, xv, x0, idx) = {
+  let noparen(expr) = {
+    if type(expr) == content and expr.func() == math.lr {
+      let children = expr.at("body").at("children")
+      children.slice(1, children.len() - 1).join()
+    } else {
+      expr
+    }
+  }
+
+  if idx == [0] or idx == 0 {
+    $fn (noparen(x0))$
+  } else if idx == [1] or idx == 1 {
+    $fn^((1)) (noparen(x0))(xv - x0)$
+  } else {
+    $frac(fn^((noparen(idx))) (noparen(x0)), idx !)(xv - x0)^noparen(idx)$
+  }
+}
+
 #let isotope(element, /*atomic mass*/a: none, /*atomic number*/z: none) = {
   $attach(upright(element), tl: #a, bl: #z)$
 }