diff --git a/README.md b/README.md
index ef100cd9..df844712 100644
--- a/README.md
+++ b/README.md
@@ -59,36 +59,41 @@ $ python setup.py install --help
 ```
 
 ## Usage
-Here is an example to solve div(alpha grad(u)) = f in a square and to plot the result
-with matplotlib (modified from ex1.cpp). Use the badge above to open this in Binder.
+This example solves the Poisson equation, $\nabla \cdot (\alpha \nabla u) = f$, in a square and plots the result
+with matplotlib (modified from `ex1.cpp`).
+Use the badge above to open this in Binder.
+
 ```python
 import mfem.ser as mfem
 
-# create sample mesh for square shape
+# Create a square mesh
 mesh = mfem.Mesh(10, 10, "TRIANGLE")
 
-# create finite element function space
+# Define the finite element function space
 fec = mfem.H1_FECollection(1, mesh.Dimension())   # H1 order=1
 fespace = mfem.FiniteElementSpace(mesh, fec)
 
-#
+# Define the essential dofs
 ess_tdof_list = mfem.intArray()
 ess_bdr = mfem.intArray([1]*mesh.bdr_attributes.Size())
 fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list)
 
-# constant coefficient (diffusion coefficient and RHS)
+# Define constants for alpha (diffusion coefficient) and f (RHS)
 alpha = mfem.ConstantCoefficient(1.0)
 rhs = mfem.ConstantCoefficient(1.0)
 
-# Note:
-#    Diffusion coefficient can be variable. To use numba-JIT compiled
-#    functio. Use the following, where x is numpy-like array.
-# @mfem.jit.scalar
-# def alpha(x):
-#     return x+1.0
-#
+"""
+Note
+-----
+In order to represent a variable diffusion coefficient, you
+must use a numba-JIT compiled function. For example:
+
+>>> @mfem.jit.scalar
+>>> def alpha(x):
+>>>     return x+1.0
+"""
 
-# define Bilinear and Linear operator
+# Define the bilinear and linear operators
 a = mfem.BilinearForm(fespace)
 a.AddDomainIntegrator(mfem.DiffusionIntegrator(alpha))
 a.Assemble()
@@ -96,38 +101,37 @@ b = mfem.LinearForm(fespace)
 b.AddDomainIntegrator(mfem.DomainLFIntegrator(rhs))
 b.Assemble()
 
-# create gridfunction, which is where the solution vector is stored
-x = mfem.GridFunction(fespace);
+# Initialize a gridfunction to store the solution vector
+x = mfem.GridFunction(fespace)
 x.Assign(0.0)
 
-# form linear equation (AX=B)
+# Form the linear system of equations (AX=B)
 A = mfem.OperatorPtr()
 B = mfem.Vector()
 X = mfem.Vector()
-a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
+a.FormLinearSystem(ess_tdof_list, x, b, A, X, B)
 print("Size of linear system: " + str(A.Height()))
 
-# solve it using PCG solver and store the solution to x
+# Solve the linear system using PCG and store the solution in x
 AA = mfem.OperatorHandle2SparseMatrix(A)
 M = mfem.GSSmoother(AA)
 mfem.PCG(AA, M, B, X, 1, 200, 1e-12, 0.0)
 a.RecoverFEMSolution(X, b, x)
 
-# extract vertices and solution as numpy array
+# Extract vertices and solution as numpy arrays
 verts = mesh.GetVertexArray()
 sol = x.GetDataArray()
 
-# plot solution using Matplotlib
-
+# Plot the solution using matplotlib
 import matplotlib.pyplot as plt
 import matplotlib.tri as tri
 
 triang = tri.Triangulation(verts[:,0], verts[:,1])
 
-fig1, ax1 = plt.subplots()
-ax1.set_aspect('equal')
-tpc = ax1.tripcolor(triang, sol, shading='gouraud')
-fig1.colorbar(tpc)
+fig, ax = plt.subplots()
+ax.set_aspect('equal')
+tpc = ax.tripcolor(triang, sol, shading='gouraud')
+fig.colorbar(tpc)
 plt.show()
 ```
 ![](https://raw.githubusercontent.com/mfem/PyMFEM/master/docs/example_image.png)