diff --git a/book/src/functions.md b/book/src/functions.md
index 7b0a50c..f1926e8 100644
--- a/book/src/functions.md
+++ b/book/src/functions.md
@@ -1,4 +1,6 @@
-# Pre-defined functions
+# Pre-defined functions & variables
+
+## Pre-defined functions
 
 The DiffSL supports the following mathematical functions that can be used in an expression:
 
@@ -7,7 +9,13 @@ The DiffSL supports the following mathematical functions that can be used in an
 * `cos(x)` - cosine of x
 * `tan(x)` - tangent of x
 * `exp(x)` - exponential of x
+* `tanh(x)` - hyperbolic tangent of x
+* `sinh(x)` - hyperbolic sine of x
+* `cosh(x)` - hyperbolic cosine of x
+* `arcsinh(x)` - inverse hyperbolic sine of x
+* `arccosh(x)` - inverse hyperbolic cosine of x
 * `log(x)` - natural logarithm of x
+* `log10(x)` - base-10 logarithm of x
 * `sqrt(x)` - square root of x
 * `abs(x)` - absolute value of x
 * `sigmoid(x)` - sigmoid function of x
@@ -15,17 +23,21 @@ The DiffSL supports the following mathematical functions that can be used in an
 
 You can use these functions as part of an expression in the DSL. For example, to define a variable `a` that is the sine of another variable `b`, you can write:
 
-```
+```diffsl
 b { 1.0 }
 a { sin(b) }
 ```
 
-# Pre-defined variables
+## Pre-defined variables
+
+There are two predefined variables in DiffSL. The first is the scalar `t` which is the current time, this allows the equations to be written as functions of time. For example
+
+```diffsl
+F_i { t + sin(t) }
+```
 
-The only predefined variable is the scalar `t` which is the current time, this allows the equations to be written as functions of time. For example
+The second is the scalar `N` which is the current model index, this allows us to write hybrid models with multiple ODE systems in the same file. For example
 
+```diffsl
+F_i { g_i[N] * x }
 ```
-F_i {
-  t + sin(t)
-}
-```
\ No newline at end of file
diff --git a/book/src/inputs_outputs.md b/book/src/inputs_outputs.md
index 9077151..1cb2b86 100644
--- a/book/src/inputs_outputs.md
+++ b/book/src/inputs_outputs.md
@@ -1,6 +1,6 @@
 # Inputs & Outputs
 
-Often it is useful to parameterize the system of equations using a set of input parameters. It is also useful to be able to extract certain variables from the system for further analysis. 
+Often it is useful to parameterize the system of equations using a set of input parameters. It is also useful to be able to extract certain outputs from the system for further analysis.
 In this section we will show how to specify inputs and outputs in the DiffSL language.
 
 ## Specifying inputs
@@ -9,18 +9,18 @@ We can override the values of any scalar variables by specifying them as input
 variables.  To do this, we add a line at the top of the code to specify that
 these are input variables:
 
-```
+```diffsl
 in { k = 1.0 }
 u { 0.1 }
 F { k * u }
 ```
 
-Here we have specified a single input parameter `k` that is used in the RHS function `F`. 
-The value of `k` is set to `1.0` in the code, but this value is only a default, and can be overridden by passing in a value at solve time.
+Here we have specified a single input parameter `k` that is used in the RHS function `F`.
+The value of `k` is set to `1.0` in the code, but this value is only there to specify the shape of the element, and will be overridden by passing in a value at solve time.
 
 We can use input parameters anywhere in the code, including in the definition of other input parameters.
 
-```
+```diffsl
 in { k = 1.0 }
 g { 2 * k }
 F { g * u }
@@ -28,24 +28,21 @@ F { g * u }
 
 or in the intial conditions of the state variables:
 
-```
+```diffsl
 in { k = 1.0 }
-u_i {
-  x = k,
-}
+u_i { x = k }
 F { u }
 ```
 
-
 ## Specifying outputs
 
 We can also specify the outputs of the system. These might be the state
 variables themselves, or they might be other variables that are calculated from
-the state variables. 
+the state variables.
 
 Here is an example where we simply output the elements of the state vector:
 
-```
+```diffsl
 u_i {
   x = 1.0,
   y = 2.0,
@@ -56,11 +53,11 @@ out_i { x, y, z }
 
 or we can derive additional outputs from the state variables:
 
-```
+```diffsl
 u_i {
   x = 1.0,
   y = 2.0,
   z = 3.0,
 }
 out { x + y + z }
-```
\ No newline at end of file
+```
diff --git a/book/src/introduction.md b/book/src/introduction.md
index 1aadcd5..48fefd6 100644
--- a/book/src/introduction.md
+++ b/book/src/introduction.md
@@ -5,10 +5,10 @@ ODE systems and is based on the idea that a ODE system can be specified by a set
 of equations of the form:
 
 $$
-M(t) \frac{\mathrm{d}\mathbf{u}}{\mathrm{d}t} = F(\mathbf{u}, t)
+M(t) \frac{\mathrm{d}\mathbf{u}}{\mathrm{d}t} = F(\mathbf{u}, \mathbf{p}, t)
 $$
 
-where \\( \mathbf{u} \\) is the vector of state variables and \\( t \\) is the time. The DSL
+where \\( \mathbf{u} \\) is the vector of state variables, \\( \mathbf{p} \\) is the vector of parameters, and \\( t \\) is the time. The DSL
 allows the user to specify the state vector \\( \mathbf{u} \\) and the RHS function \\( F \\).
 Optionally, the user can also define the derivative of the state vector \\( \mathrm{d}\mathbf{u}/\mathrm{d}t \\)
 and the mass matrix \\( M \\) as a function of \\( \mathrm{d}\mathbf{u}/\mathrm{d}t \\) (note that this function should be linear!).
@@ -18,7 +18,7 @@ scalars and vectors of the users that are required to calculate \\( F \\) and \\
 
 ## A Simple Example
 
-To illustrate the language, consider the following simple example of a logistic growth model:
+To illustrate the language, consider the following example of a logistic growth model:
 
 $$
 \frac{\mathrm{d}N}{\mathrm{d}t} = r N (1 - N/K)
@@ -28,7 +28,7 @@ where \\( N \\) is the population, \\( r \\) is the growth rate, and \\( k \\) i
 
 To specify this model in DiffSL, we can write:
 
-```
+```diffsl
 in_i { r = 1.0, k = 10.0 }
 u_i {
   N = 0.0
@@ -43,4 +43,4 @@ out_i {
 
 Here, we define the input parameters for our model as a vector `in` with the growth rate `r` and the carrying capacity `k`. We then define the state vector `u_i` with the population `N` initialized to `0.0`. Next, we define the RHS function `F_i` as the logistic growth equation. Finally, we define the output vector `out_i` with the population `N`.
 
-*For an LLM-oriented summary and index of this documentation, see [llms.txt](./llms.txt).*
\ No newline at end of file
+*For an LLM-oriented summary and index of this documentation, see [llms.txt](./llms.txt).*
diff --git a/book/src/odes.md b/book/src/odes.md
index 0f359a6..663def1 100644
--- a/book/src/odes.md
+++ b/book/src/odes.md
@@ -12,21 +12,20 @@ $$
 where \\( \mathbf{u} \\) is the vector of state variables, \\( \mathbf{u}_0 \\) is the initial condition, \\( \mathbf{p} \\) is the parameter vector, \\( F_N \\) is the model-indexed RHS function, and \\( M \\) is the mass matrix.
 The DSL allows the user to specify the state vector \\( \mathbf{u} \\), parameter vector \\( \mathbf{p} \\), and RHS function \\( F_N \\).
 
-Optionally, the user can also define the derivative of the state vector \\( \mathrm{d}\mathbf{u}/\mathrm{d}t \\) and the mass matrix \\( M \\) as a function of \\( \mathrm{d}\mathbf{u}/\mathrm{d}t \\) 
-(note that this function should be linear!). 
-The user is also free to define an arbitrary number of intermediate tensors that are required to calculate \\( F \\) and \\( M \\). 
-
+Optionally, the user can also define the derivative of the state vector \\( \mathrm{d}\mathbf{u}/\mathrm{d}t \\) and the mass matrix \\( M \\) as a function of \\( \mathrm{d}\mathbf{u}/\mathrm{d}t \\)
+(note that this function should be linear!).
+The user is also free to define an arbitrary number of intermediate tensors that are required to calculate \\( F \\) and \\( M \\).
 
 ## Defining the state vector
 
 To define the state vector \\(\mathbf{u} \\), we create a special
 vector tensor `u_i`. Note the particular name `u` is used to indicate that this
-is the state vector.
+is the state vector, and cannot be used for any othr purpose.
 
 The values that we use for `u_i` are the initial values of the state variables
 at \\( t=0 \\), so an initial condition \\( \mathbf{u}(t=0) = [x(0), y(0), z(0)] = [1, 0, 0] \\) is defined as:
 
-```
+```diffsl
 u_i {
   x = 1,
   y = 0,
@@ -34,18 +33,18 @@ u_i {
 }
 ```
 
-Since we will often use the individual elements of the state vector in the RHS function, it is useful to define them as separate variables as well.
+Since we will often use the individual elements of the state vector in the RHS function, it is useful to define them using labels (`x`, `y`, `z`) so we can use them later.
 
 The state tensor `u` must be either be a scalar or a vector. So, if we only had a single state variable, we could define it as a scalar without the index `i` as follows:
 
-```
+```diffsl
 u { x = 1 }
 ```
 
 We can optionally define the time derivatives of the state variables,
 \\( \mathbf{\dot{u}} \\) as well:
 
-```
+```diffsl
 dudt_i {
   dxdt = 1,
   dydt = 0,
@@ -53,18 +52,15 @@ dudt_i {
 }
 ```
 
-Here the initial values of the time derivatives are given.
-In many cases any values can be given here as the time derivatives of the state variables are calculated from the RHS.
-However, if there are any algebraic variables in the equations then these values can be used 
-used as a starting point to calculate a set of consistent initial values for the
-state variables.
+Here the initial values of the time derivatives are given. These values are simply placeholders and are only used
+to specify the shape of each element, for example they do not need to be consistent with the RHS function `F_i` that we will define later.
 
 Note that there is no need to define `dudt` if you do not define a mass matrix \\( M \\).
 
 ## Defining the ODE system equations
 
 We now define the right-hand-side  function \\( F \\) that we want to solve, using the
-variables that we have defined earlier. We do this by defining a vector variable
+variables that we have defined earlier. We do this by defining a vector (or scalar)
 `F_i` that corresponds to the RHS of the equations.
 
 For example, to define a simple system of ODEs:
@@ -80,25 +76,25 @@ $$
 
 We write:
 
-```
+```diffsl
 u_i { x = 1, y = 0 }
 F_i { y, -x }
 ```
 
-## Using the Ode system index `N`
+## The system index `N`
 
-**Note: hybrid ODEs using `N` and `reset` are not yet supported in the current version of diffsol or pydiffsol, please check back for updates**
+**Note: hybrid ODEs using `N` and `reset` are only supported in the current version of diffsol or pydiffsol, please check back for updates**
 
 The index `N` is used to define multiple ODE systems in the same file, indexed by the non-negative integer `N`.
-This allows us to define multiple ODE systems in the same file. Combined with the stop and reset functions (see below), 
+This allows us to define multiple ODE systems in the same file. Combined with the stop and reset functions (see below),
 this also allows us to define hybrid switching systems where we can switch between different ODE systems at different times during the simulation.
 
-The model index `N` can be used in any of the equations, and it can be used to index into any of the tensors that we have defined.
+The model index `N` can be used in any of the equations, and it can be used to index into any vectors that we have defined.
 
 For example, we can define two ODE systems, one with exponential growth and one with exponential decay, by
 defining a `g_i` vector variable that contains the growth/decay rate for each system, and then using this variable in the definition of the RHS function `F_i`:
 
-```
+```diffsl
 u_i { x = 1 }
 g_i { 0.1, -0.1 }
 F_i { g_i[N] * x }
@@ -109,44 +105,25 @@ F_i { g_i[N] * x }
 The `stop` and `reset` functions are standard tensors that can be used to stop and reset the ODE system at a given time, as long as the runtime supports this feature.
 
 The `reset` tensor should be the same shape (i.e. a vector with the same number of elements) as the state vector `u_i`, whereas the `stop` tensor can be any length vector.
-During the solve, whenever any element of the `stop` tensor is equal to zero, the ODE system will stop and, if the runtime supports hybrid models, the state of the 
+During the solve, whenever any element of the `stop` tensor is equal to zero, the ODE system will stop and, if the runtime supports hybrid models, the state of the
 ODE system will be reset to the values defined in the `reset` tensor, and the ODE system will continue to solve from there.
 
 The classic example of this type of hybrid model is a bouncing ball, where the ODE system is stopped when the ball hits the ground (\(x \leq 0\)), and then the state of the system is reset to a new value that corresponds to the ball bouncing back up.
 
-```
-u_i {
- x = 0,
- v = 10,
-}
-F_i {
- v,
- -9.81,
-}
-stop {
- x,
-}
-reset {
- 0,
- -0.8 * v,
-}
+```diffsl
+u_i { x = 0, v = 10 }
+F_i { v, -9.81 }
+stop { x }
+reset { 0, -0.8 * v }
 ```
 
 Another example is a system that is periodically forced, for example dosing into a pharmacokinetic model, where the ODE system is stopped at regular intervals (e.g. every 24 hours), and the state of the system is reset to a new value that corresponds to the dose being administered.
 
-```
-u_i {
- x = 0,
-}
-F_i {
- -0.1 * x,
-}
-stop {
- t - 24,
-}
-reset {
- x + 10,
-}
+```diffsl
+u { x = 0 }
+F { -0.1 * x }
+stop { t - 24 }
+reset { x + 10 }
 ```
 
 ## Hybrid models with multiple ODE systems
@@ -155,18 +132,17 @@ The `stop` and `reset` functions can also be used to switch between different OD
 
 For example, we can define a system that starts with exponential growth, and then switches to exponential decay after 24 hours, and then switches back to exponential growth after another 24 hours, by defining two ODE systems as before, and then using the `stop` and `reset` functions to switch between them:
 
-
-```
-u_i { x = 1 }
+```diffsl
+u { x = 1 }
 g_i { 0.1, -0.1 }
-F_i { g_i[N] * x }
+F { g_i[N] * x }
 stop { t - 48, t - 24 }
 reset { 1 }
 ```
 
 ## Defining the mass matrix
 
-We can also define a mass matrix \\( M \\) by defining a vector variable `M_i` which is the product of the mass matrix with the time derivative of the state vector \\( M \mathbf{\dot{u}} \\). 
+We can also define a mass matrix \\( M \\) by defining a vector variable `M_i` which is the product of the mass matrix with the time derivative of the state vector \\( M \mathbf{\dot{u}} \\).
 This is optional, and if not defined, the mass matrix is assumed to be the identity
 matrix.
 
@@ -186,21 +162,9 @@ $$
 
 We write:
 
+```diffsl
+u_i { x = 1, y = 0 }
+dudt_i { dxdt = 0, dydt = 1 }
+M_i { dxdt, 0 }
+F_i { x, y-x }
 ```
-u_i {
- x = 1,
- y = 0,
-}
-dudt_i {
- dxdt = 0,
- dydt = 1,
-}
-M_i {
- dxdt,
- 0,
-}
-F_i {
- x,
- y-x,
-}
-```
\ No newline at end of file
diff --git a/book/src/operations.md b/book/src/operations.md
index ff38274..716dce8 100644
--- a/book/src/operations.md
+++ b/book/src/operations.md
@@ -1,13 +1,13 @@
 # Tensor Operations
 
-We can use standard algebraic operations on variables. To refer to previously
-defined variables, we use the variable name, making sure to use the correct
+We can use standard algebraic operations on tensors. To refer to previously
+defined variables, we use the tensor name, making sure to use the correct
 subscript if it is a vector or tensor.
 
-For example, to define a scalar variable \\( a \\) as the sum of two other scalar
-variables \\( b \\) and \\( c \\), we write:
+For example, to define a scalar \\( a \\) as the sum of two other scalars
+\\( b \\) and \\( c \\), we write:
 
-```
+```diffsl
 b { 1.0 }
 c { 2.0 }
 a { b + c }
@@ -15,10 +15,10 @@ a { b + c }
 
 The scalar `a` will therefore be equal to 3.0.
 
-To define a vector variable \\( \mathbf{v} \\) as the sum of two other vector
-variables \\( \mathbf{u} \\) and \\( \mathbf{w} \\), we write:
+To define a vector \\( \mathbf{v} \\) as the sum of two other vectors
+\\( \mathbf{u} \\) and \\( \mathbf{w} \\), we write:
 
-```
+```diffsl
 u_i { 1.0, 2.0 }
 w_i { 3.0, 4.0 }
 v_i { u_i + w_i }
@@ -27,7 +27,7 @@ v_i { u_i + w_i }
 Notice that the index of the vectors within the expression must match the index of the output vector `v_i`.
 So if we defined `v_a` instead of `v_i`, the expression would be:
 
-```
+```diffsl
 v_a { u_a + w_a }
 ```
 
@@ -38,7 +38,7 @@ For higher-dimensional tensors, the order of the indices in the expression relat
 For example, we can use the indices to define a translation of a matrix. Here we define a new matrix \\( C \\) that is the sum of \\( A \\) and \\( B^T \\),
 where \\( B^T \\) is the transpose of \\( B \\)
 
-```
+```diffsl
 C_ij { A_ij + B_ji }
 ```
 
@@ -48,7 +48,7 @@ Notice that the indices of \\( B^T \\) are reversed in the expression compared t
 
 Broadcasting is supported in the language, so you can perform element-wise operations on tensors of different shapes. For example, the following will define a new vector \\( \mathbf{d} \\) that is the sum of \\( \mathbf{a} \\) and a scalar \\( k \\):
 
-```
+```diffsl
 a_i { 1.0, 2.0 }
 k { 3.0 }
 d_i { a_i + k }
@@ -64,17 +64,19 @@ DiffSL broadcasting is index-based:
 
 This means broadcasting behavior depends on the index labels used in the expression, not only on raw tensor shapes.
 
-For example, the following is valid because `b` is indexed by `j`, so it broadcasts across the `i` axis:
+For example, the following example defined a 3x2 matrix `A_ij` and a vector `b_i` of length 3. Given these shapes, the following definition
+of `c_ij` is valid because `b` is indexed by `j` in the expression `A_ij + b_j`. Since the number of *rows* in `b` matches the *columns* of `A`,
+this is a valid broadcasting operation.
 
-```
+```diffsl
 A_ij { (0:3, 0:2): 1.0 }
 b_i { (0:2): 1.0 }
 c_ij { A_ij + b_j }
 ```
 
-The following is invalid because `b_i` is aligned to the first axis (`i`) and cannot be broadcast to the matrix shape:
+However, if `b` is instead indexed by `i`, the broadcasting is invalid because `b_i` is aligned to the *rows* (the `i` axis) of `A`, which are a different length, so the broadcasting cannot be performed and a compiler error will be raised.
 
-```
+```diffsl
 A_ij { (0:3, 0:2): 1.0 }
 b_i { (0:2): 1.0 }
 c_ij { A_ij + b_i }
@@ -82,52 +84,52 @@ c_ij { A_ij + b_i }
 
 ## Contractions
 
-The DiffSL indexing notation allows for tensor contractions, which are operations that sum over one or more indices. 
-The rule used is that any indices that do not appear in the output tensor will be summed over. 
+The DiffSL indexing notation allows for tensor contractions, which are operations that sum over one or more indices.
+The rule used is that any indices that do not appear in the output tensor will be summed over.
 
-For example, the following defines a new vector \\( \mathbf{v} \\) that is the sum of the rows of a matrix \\( A \\):
+For example, the following defines a new vector \\( \mathbf{v} \\) that is the sum of each individual row of matrix \\( A \\) (i.e. the sum is taken over the `j` index of `A`):
 
-```
+```diffsl
 v_i { A_ij }
 ```
 
 The above expression sums over the `j` index of the matrix `A`, resulting in a vector `v` where each element `v_i` is the sum of the elements in the `i`-th row of `A`.
-At the moment only 2d to 1d contractions are supported in order to enable matrix-vector multiplication, please open an issue if you need more general contraction support.
+At the moment only 2d to 1d contractions are supported in order to enable matrix-vector multiplication, please comment on [this issue](https://github.com/martinjrobins/diffsl/issues/76) if you need more general contraction support.
 
-We can also define a matrix-vector multiplication, the following will define a new vector \\( \mathbf{v} \\) that is
+Using a contraction we can define the popular matrix-vector multiplication operation. The following will define a new vector \\( \mathbf{v} \\) that is
 the result of a matrix-vector multiplication of a matrix \\( A \\) and a vector \\( \mathbf{u} \\):
 
-```
+```diffsl
 v_i { A_ij * u_j }
 ```
 
-This operation is actually a combination of a broadcast `A_ij * u_j`, followed by a contraction over the `j` index, 
-the `A_ij * u_j` expression broadcasts the vector `u` to the same shape as `A`, forming a new 2D tensor, and 
-the output vector `v_i` implicitly sums over the missing `j` index to form the final output vector. 
+This operation is actually a combination of a broadcast `A_ij * u_j`, followed by a contraction over the `j` index,
+the `A_ij * u_j` expression broadcasts the vector `u` to the same shape as `A`, forming a new 2D tensor, and
+the output vector `v_i` implicitly sums over the missing `j` index to form the final output vector.
 To illustrate this further, lets manually break this matrix-vector multiplication into two steps using an intermediary tensor `M_ij`:
 
-```
+```diffsl
 M_ij { A_ij * u_j }
 v_i { M_ij }
 ```
 
 The first step calculates the element-wise product of `A` and `u` using broadcasting into the 2D tensor `M`, and the second step uses a contraction to sum over the `j` index to form the output vector `v`.
 
-
 ## Indexing
 
-Indexing a 1D dense tensor (vector) is supported using square brackets. You can use either single indexing to extract a single element, or range indexing to extract a sub-vector.
+Indexing a 1D dense tensor (vector) is supported using square brackets. If you need more general indexing capabilities, please comment on [this issue](https://github.com/martinjrobins/diffsl/issues/77).
 
+You can use either single indexing to extract a single element, or range indexing to extract a sub-vector.
 For example, to extract the third element of a vector \\( \mathbf{a} \\) and assign it to a scalar \\( r \\), you can write:
 
-```
+```diffsl
 a_i { 0.0, 1.0, 2.0, 3.0 }
 r { a_i[2] }
 ```
 
 To extract a sub-vector containing the second and third elements of \\( \mathbf{a} \\) and assign it to a new vector \\( \mathbf{r} \\), you can write:
 
-```
+```diffsl
 a_i { 0.0, 1.0, 2.0, 3.0 }
 r_i { a_i[1:3] }
 ```
diff --git a/book/src/tensors.md b/book/src/tensors.md
index 2c99fb1..fe4fa7b 100644
--- a/book/src/tensors.md
+++ b/book/src/tensors.md
@@ -5,7 +5,7 @@ These tensors can be dense, sparse or diagonal, and the compiler will try to cho
 
 The simplest tensor is a 0th dimensional scalar. For example, to define a scalar variable \\( k \\) with value 1.0, we write:
 
-```
+```diffsl
 k { 1.0 }
 ```
 
@@ -14,7 +14,7 @@ and here `k` only has a single element, which is a constant value 1.
 
 Lets now define a 1-dimensional vector variable \\( \mathbf{v} \\) with 3 elements:
 
-```
+```diffsl
 v_i {
  1.0,
  2.0,
@@ -25,39 +25,50 @@ v_i {
 The list of elements within a tensor are deliniated with a comma `,` and the trailing comma at the end of the list is optional.
 Whitespace is ignored so you can also write this tensor all on the one line:
 
-```
+```diffsl
 v_i { 1.0, 2.0, 3.0 }
 ```
 
 ## Subscripts
 
-In the previous vector `v_i`, the subscript `_i` indicates that this is a 1D vector. Each subscript is a single character, 
+In the previous vector `v_i`, the subscript `_i` indicates that this is a 1D vector. Each subscript is a single character,
 and the number of subscripts indicates the number of dimensions of the tensor. You can use any character for the subscript,
 
-```
+```diffsl
 v_x { 1.0, 2.0, 3.0 }
 w_a { 2.0, 3.0 }
 v_i { 1.0, 2.0, 3.0, 4.0 }
 ```
 
+If we wanted to define a 2D matrix, we would need to use two subscripts, for example:
+
+```diffsl
+A_ij { (0:2, 0:3): 1.0 }
+```
+
+If we wanted to define a 0D scalar, we would simply not use any subscripts:
+
+```diffsl
+k { 1.0 }
+```
+
 ## Ranges
 
 Each element of a tensor can optionally give a *range* of index numbers, which is used by the compiler to determine the extent of each element.
 This is useful when defining higher dimensional tensors, such as matrices. For example, to define a 2x3 matrix \\( A \\) with all elements set to `1.0`, we write:
 
-```
+```diffsl
 A_ij {
  (0:2, 0:3): 1.0,
 }
 ```
 
-
 Note the two subscript to indicate that this is a 2D tensor. The size of the
 single element is given in the brackets, we have two ranges `0:2` and `0:3` that correspond to the two dimensions of the matrix.
 
 Here is another example of a 4x2 matrix \\( B \\) with rows 0 to 2 set to `1.0` and rows 3 to 4 set to `2.0`:
 
-```
+```diffsl
 A_ij {
  (0:2, 0:3): 1.0,
  (3:4, 0:3): 2.0,
@@ -66,7 +77,7 @@ A_ij {
 
 For specifying a single index, you can simply write the index number without the colon, for example to define a 3x3 identity matrix $I$:
 
-```
+```diffsl
 I_ij {
  (0, 0): 1.0,
  (1, 1): 1.0,
@@ -81,7 +92,7 @@ Notice also that we have not defined all the elements of the matrix, only the no
 
 Finally, you can also use the `..` operator to specify a *diagonal* range of indices. For example, to define a 3x3 matrix \\( D \\) with the diagonal elements set to `1.0`:
 
-```
+```diffsl
 D_ij {
  (0..2, 0..2): 1.0,
 }
@@ -91,7 +102,7 @@ D_ij {
 
 We can automatically define a sparse matrix \\( B \\) by simply specifying the non-zero elements:
 
-```
+```diffsl
 B_ij {
  (0, 0): 1.0,
  (0, 1): 2.0,
@@ -100,20 +111,23 @@ B_ij {
 ```
 
 The compiler will infer that this is a 2x2 matrix, and will automatically represent it as a sparse matrix.
-We can force the compiler to use a dense representation by specifying the zeros explicitly:
 
-```
-B_ij {
- (0, 0): 1.0,
- (0, 1): 2.0,
- (1, 0): 0.0,
- (1, 1): 3.0,
+Any zero-value elements will automatically be filtered from the matrix, this is useful for specifying the sizes of
+sparse matrices if the last row or column is zero. For example, to define a 3x3 matrix \\( C \\) with the first two rows set to `1.0`:
+
+```diffsl
+C_ij {
+ (0:2, 0:3): 1.0,
+ (2, 2): 0.0,
 }
 ```
 
+Without the zero element, the compiler would infer that this is a 2x3 matrix, but with the zero element it will infer that this is a 3x3 matrix with the last row set to zero.
+The sparsity of the matrix will be preserved in the internal representation, so the zero elements will not take up any memory or computational resources.
+
 As well as specifying a sparse matrix, we can also define a diagonal matrix by specifying the diagonal elements:
 
-```
+```diffsl
 D_ij {
  (0, 0): 1.0,
  (1, 1): 2.0,
@@ -123,13 +137,13 @@ D_ij {
 
 The compiler will infer that this is a 3x3 matrix, and will automatically represent it as a diagonal matrix.
 
-## Labels
+## Element Labels
 
 Each element of a tensor can optionally be given a name or *label* that can be used to refer to the element in expressions.
 
 For example, to define a vector with two named elements:
 
-```
+```diffsl
 v_i {
  x = 1.0,
  y = 2.0,
@@ -138,9 +152,7 @@ v_i {
 
 Here we have defined a single tensor `v_i` with two named elements `x` and `y`. We can then refer to the individual elements in expressions, where they will act as if they were separate variables:
 
-```
+```diffsl
 v_i { x = 1.0, y = 2.0 }
 w_i { 2 * y, 3 * x }
 ```
-
-