diff --git a/demo/sofa.py b/demo/sofa.py
new file mode 100644
index 00000000..00a7c619
--- /dev/null
+++ b/demo/sofa.py
@@ -0,0 +1,63 @@
+'''
+This script calculates the constant in Gerver's solution to the moving sofa
+problem.
+
+See Finch, S. R. "Moving Sofa Constant." §8.12 in Mathematical Constants.
+Cambridge, England: Cambridge University Press, pp. 519-523, 2003.
+'''
+
+from mpmath import cos, sin, pi, quad, findroot, mp
+
+mp.prec = 113
+
+eqs = [lambda A, B, φ, θ: (A*(cos(θ) - cos(φ)) - 2*B*sin(φ)
+                           + (θ - φ - 1)*cos(θ) - sin(θ) + cos(φ) + sin(φ)),
+       lambda A, B, φ, θ: (A*(3*sin(θ) + sin(φ)) - 2*B*cos(φ)
+                           + 3*(θ - φ - 1)*sin(θ) + 3*cos(θ) - sin(φ) + cos(φ)),
+       lambda A, B, φ, θ: A*cos(φ) - (sin(φ) + 0.5 - 0.5*cos(φ) + B*sin(φ)),
+       lambda A, B, φ, θ: ((A + pi/2 - φ - θ) - (B - (θ - φ)*(1 + A)/2
+                                                 - 0.25*(θ - φ)**2))]
+A, B, φ, θ = findroot(eqs, (0, 0, 0, 0))
+
+def r(α):
+    if 0 <= α < φ:
+        return 0.5
+    if φ <= α < θ:
+        return (1 + A + α - φ)/2
+    if θ <= α < pi/2 - θ:
+        return A + α - φ
+    return B - (pi/2 - α - φ)*(1 + A)/2 - (pi/2 - α - φ)**2/4
+
+s = lambda α: 1 - r(α)
+
+def u(α):
+    if φ <= α < θ:
+        return B - (α - φ)*(1 + A)/2 - (α - φ)**2/4
+    return A + pi/2 - φ - α
+
+def du(α):
+    if φ <= α < θ:
+        return -(1 + A)/2 - (α - φ)/2
+    return -1
+
+def y(α, f):
+    if α > pi/2 - θ:
+        i = [0, φ, θ, pi/2 - θ, α]
+    elif α > θ:
+        i = [0, φ, θ, α]
+    elif α > φ:
+        i = [0, φ, α]
+    else:
+        i = i = [0, α]
+    return 1 - quad(lambda x: f(x)*sin(x), i)
+
+y1 = lambda α: y(α, r)
+y2 = lambda α: y(α, s)
+y3 = lambda α: y2(α) - u(α)*sin(α)
+
+S1 = quad(lambda x: y1(x)*r(x)*cos(x), [0, φ, θ, pi/2 - θ, pi/2 - φ])
+S2 = quad(lambda x: y2(x)*s(x)*cos(x), [0, φ, θ])
+S3 = quad(lambda x: y3(x)*(u(x)*sin(x) - du(x)*cos(x) - s(x)*cos(x)),
+          [φ, θ, pi/4])
+
+print(2*(S1 + S2 + S3))