diff --git a/docs/source/expression_tree/state_vector.rst b/docs/source/expression_tree/state_vector.rst
index 4a0792b3b7..4370f6e980 100644
--- a/docs/source/expression_tree/state_vector.rst
+++ b/docs/source/expression_tree/state_vector.rst
@@ -3,3 +3,6 @@ State Vector
 
 .. autoclass:: pybamm.StateVector
   :members:
+
+.. autoclass:: pybamm.StateVectorDot
+  :members:
diff --git a/docs/source/expression_tree/variable.rst b/docs/source/expression_tree/variable.rst
index 0e610f5f28..510e366cf1 100644
--- a/docs/source/expression_tree/variable.rst
+++ b/docs/source/expression_tree/variable.rst
@@ -4,5 +4,9 @@ Variable
 .. autoclass:: pybamm.Variable
   :members:
 
+.. autoclass:: pybamm.VariableDot
+  :members:
+
 .. autoclass:: pybamm.ExternalVariable
   :members:
+
diff --git a/examples/notebooks/expression_tree/expression-tree.ipynb b/examples/notebooks/expression_tree/expression-tree.ipynb
index 9bcfeb3194..bf28c8fe2e 100644
--- a/examples/notebooks/expression_tree/expression-tree.ipynb
+++ b/examples/notebooks/expression_tree/expression-tree.ipynb
@@ -64,7 +64,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "We can also calculate the expression tree representing the gradient of the equation with respect to $t$ (which is of course simply the scalar value 1),"
+    "We can also calculate the expression tree representing the gradient of the equation with respect to $t$,"
    ]
   },
   {
@@ -84,7 +84,7 @@
     "![](expression_tree2.png)\n",
     "\n",
     "\n",
-    "...and evaluate this expression, which will again give 1."
+    "...and evaluate this expression,"
    ]
   },
   {
@@ -95,7 +95,7 @@
     {
      "data": {
       "text/plain": [
-       "1.0"
+       "array([[-11.]])"
       ]
      },
      "execution_count": 4,
@@ -104,7 +104,7 @@
     }
    ],
    "source": [
-    "diff_wrt_equation.evaluate(1, np.array([2]))"
+    "diff_wrt_equation.evaluate(t=1, y=np.array([2]), y_dot=np.array([2]))"
    ]
   },
   {
@@ -202,6 +202,13 @@
     "\n",
     "After the third stage, our expression tree is now able to be evaluated by one of the solver classes. Note that we have used a single equation above to illustrate the different types of expression trees in PyBaMM, but any given models will consist of many RHS or algebraic equations, along with boundary conditions. See [here](https://github.com/pybamm-team/PyBaMM/blob/master/examples/notebooks/add-model.ipynb) for more details of PyBaMM models."
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {
@@ -220,7 +227,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.3"
+   "version": "3.6.7"
   }
  },
  "nbformat": 4,
diff --git a/examples/notebooks/expression_tree/expression_tree2.png b/examples/notebooks/expression_tree/expression_tree2.png
index 949690cfc0..09f6e2900e 100644
Binary files a/examples/notebooks/expression_tree/expression_tree2.png and b/examples/notebooks/expression_tree/expression_tree2.png differ
diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index 9ca4cf1ee2..6fd55219f2 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -77,11 +77,12 @@ def version(formatted=False):
 from .expression_tree.parameter import Parameter, FunctionParameter
 from .expression_tree.broadcasts import *
 from .expression_tree.scalar import Scalar
-from .expression_tree.variable import Variable, ExternalVariable
+from .expression_tree.variable import Variable, ExternalVariable, VariableDot
+from .expression_tree.variable import VariableBase
 from .expression_tree.independent_variable import *
 from .expression_tree.independent_variable import t
 from .expression_tree.vector import Vector
-from .expression_tree.state_vector import StateVector
+from .expression_tree.state_vector import StateVectorBase, StateVector, StateVectorDot
 
 from .expression_tree.exceptions import *
 
diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index 8f4e4ea3a5..b3b371c147 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -515,6 +515,7 @@ def process_rhs_and_algebraic(self, model):
             equations) and processed_concatenated_algebraic
 
         """
+
         # Discretise right-hand sides, passing domain from variable
         processed_rhs = self.process_dict(model.rhs)
 
@@ -856,6 +857,13 @@ def _process_symbol(self, symbol):
             disc_children = [self.process_symbol(child) for child in symbol.children]
             return symbol._function_new_copy(disc_children)
 
+        elif isinstance(symbol, pybamm.VariableDot):
+            return pybamm.StateVectorDot(
+                *self.y_slices[symbol.get_variable().id],
+                domain=symbol.domain,
+                auxiliary_domains=symbol.auxiliary_domains
+            )
+
         elif isinstance(symbol, pybamm.Variable):
             # Check if variable is a standard variable or an external variable
             if any(symbol.id == var.id for var in self.external_variables.values()):
diff --git a/pybamm/expression_tree/array.py b/pybamm/expression_tree/array.py
index 348fbed616..43a8864eb7 100644
--- a/pybamm/expression_tree/array.py
+++ b/pybamm/expression_tree/array.py
@@ -98,6 +98,6 @@ def new_copy(self):
             self.entries_string,
         )
 
-    def _base_evaluate(self, t=None, y=None, u=None):
+    def _base_evaluate(self, t=None, y=None, y_dot=None, u=None):
         """ See :meth:`pybamm.Symbol._base_evaluate()`. """
         return self._entries
diff --git a/pybamm/expression_tree/binary_operators.py b/pybamm/expression_tree/binary_operators.py
index 392c1641ab..cadd01267f 100644
--- a/pybamm/expression_tree/binary_operators.py
+++ b/pybamm/expression_tree/binary_operators.py
@@ -13,7 +13,7 @@ def is_scalar_zero(expr):
     Utility function to test if an expression evaluates to a constant scalar zero
     """
     if expr.is_constant():
-        result = expr.evaluate_ignoring_errors()
+        result = expr.evaluate_ignoring_errors(t=None)
         return isinstance(result, numbers.Number) and result == 0
     else:
         return False
@@ -24,7 +24,7 @@ def is_matrix_zero(expr):
     Utility function to test if an expression evaluates to a constant matrix zero
     """
     if expr.is_constant():
-        result = expr.evaluate_ignoring_errors()
+        result = expr.evaluate_ignoring_errors(t=None)
         return (issparse(result) and result.count_nonzero() == 0) or (
             isinstance(result, np.ndarray) and np.all(result == 0)
         )
@@ -32,12 +32,12 @@ def is_matrix_zero(expr):
         return False
 
 
-def is_one(expr):
+def is_scalar_one(expr):
     """
     Utility function to test if an expression evaluates to a constant scalar one
     """
     if expr.is_constant():
-        result = expr.evaluate_ignoring_errors()
+        result = expr.evaluate_ignoring_errors(t=None)
         return isinstance(result, numbers.Number) and result == 1
     else:
         return False
@@ -162,21 +162,21 @@ def _binary_new_copy(self, left, right):
         "Default behaviour for new_copy"
         return self.__class__(left, right)
 
-    def evaluate(self, t=None, y=None, u=None, known_evals=None):
+    def evaluate(self, t=None, y=None, y_dot=None, u=None, known_evals=None):
         """ See :meth:`pybamm.Symbol.evaluate()`. """
         if known_evals is not None:
             id = self.id
             try:
                 return known_evals[id], known_evals
             except KeyError:
-                left, known_evals = self.left.evaluate(t, y, u, known_evals)
-                right, known_evals = self.right.evaluate(t, y, u, known_evals)
+                left, known_evals = self.left.evaluate(t, y, y_dot, u, known_evals)
+                right, known_evals = self.right.evaluate(t, y, y_dot, u, known_evals)
                 value = self._binary_evaluate(left, right)
                 known_evals[id] = value
                 return value, known_evals
         else:
-            left = self.left.evaluate(t, y, u)
-            right = self.right.evaluate(t, y, u)
+            left = self.left.evaluate(t, y, y_dot, u)
+            right = self.right.evaluate(t, y, y_dot, u)
             return self._binary_evaluate(left, right)
 
     def _evaluate_for_shape(self):
@@ -252,8 +252,12 @@ def _binary_simplify(self, left, right):
         if is_scalar_zero(right):
             return pybamm.Scalar(1)
 
-        # anything to the power of one is itself
+        # zero to the power of anything is zero
         if is_scalar_zero(left):
+            return pybamm.Scalar(0)
+
+        # anything to the power of one is itself
+        if is_scalar_one(right):
             return left
 
         return self.__class__(left, right)
@@ -425,9 +429,9 @@ def _binary_simplify(self, left, right):
             return zeros_of_shape(shape)
 
         # anything multiplied by a scalar one returns itself
-        if is_one(left):
+        if is_scalar_one(left):
             return right
-        if is_one(right):
+        if is_scalar_one(right):
             return left
 
         return pybamm.simplify_multiplication_division(self.__class__, left, right)
@@ -549,7 +553,7 @@ def _binary_simplify(self, left, right):
                 return pybamm.Array(np.inf * np.ones(left.shape_for_testing))
 
         # anything divided by one is itself
-        if is_one(right):
+        if is_scalar_one(right):
             return left
 
         return pybamm.simplify_multiplication_division(self.__class__, left, right)
@@ -622,9 +626,9 @@ def _binary_simplify(self, left, right):
             return zeros_of_shape(shape)
 
         # anything multiplied by a scalar one returns itself
-        if is_one(left):
+        if is_scalar_one(left):
             return right
-        if is_one(right):
+        if is_scalar_one(right):
             return left
 
         return pybamm.simplify_multiplication_division(self.__class__, left, right)
diff --git a/pybamm/expression_tree/concatenations.py b/pybamm/expression_tree/concatenations.py
index e1225bf03a..4a74de84cf 100644
--- a/pybamm/expression_tree/concatenations.py
+++ b/pybamm/expression_tree/concatenations.py
@@ -54,7 +54,7 @@ def _concatenation_evaluate(self, children_eval):
         else:
             return self.concatenation_function(children_eval)
 
-    def evaluate(self, t=None, y=None, u=None, known_evals=None):
+    def evaluate(self, t=None, y=None, y_dot=None, u=None, known_evals=None):
         """ See :meth:`pybamm.Symbol.evaluate()`. """
         children = self.cached_children
         if known_evals is not None:
@@ -62,14 +62,14 @@ def evaluate(self, t=None, y=None, u=None, known_evals=None):
                 children_eval = [None] * len(children)
                 for idx, child in enumerate(children):
                     children_eval[idx], known_evals = child.evaluate(
-                        t, y, u, known_evals
+                        t, y, y_dot, u, known_evals
                     )
                 known_evals[self.id] = self._concatenation_evaluate(children_eval)
             return known_evals[self.id], known_evals
         else:
             children_eval = [None] * len(children)
             for idx, child in enumerate(children):
-                children_eval[idx] = child.evaluate(t, y, u)
+                children_eval[idx] = child.evaluate(t, y, y_dot, u)
             return self._concatenation_evaluate(children_eval)
 
     def new_copy(self):
diff --git a/pybamm/expression_tree/functions.py b/pybamm/expression_tree/functions.py
index cd6b50bf01..f6f87b5c78 100644
--- a/pybamm/expression_tree/functions.py
+++ b/pybamm/expression_tree/functions.py
@@ -152,19 +152,20 @@ def _function_jac(self, children_jacs):
 
         return jacobian
 
-    def evaluate(self, t=None, y=None, u=None, known_evals=None):
+    def evaluate(self, t=None, y=None, y_dot=None, u=None, known_evals=None):
         """ See :meth:`pybamm.Symbol.evaluate()`. """
         if known_evals is not None:
             if self.id not in known_evals:
                 evaluated_children = [None] * len(self.children)
                 for i, child in enumerate(self.children):
                     evaluated_children[i], known_evals = child.evaluate(
-                        t, y, u, known_evals=known_evals
+                        t, y, y_dot, u, known_evals=known_evals
                     )
                 known_evals[self.id] = self._function_evaluate(evaluated_children)
             return known_evals[self.id], known_evals
         else:
-            evaluated_children = [child.evaluate(t, y, u) for child in self.children]
+            evaluated_children = [child.evaluate(t, y, y_dot, u)
+                                  for child in self.children]
             return self._function_evaluate(evaluated_children)
 
     def _evaluate_for_shape(self):
diff --git a/pybamm/expression_tree/independent_variable.py b/pybamm/expression_tree/independent_variable.py
index 8c11288d98..b0d8990931 100644
--- a/pybamm/expression_tree/independent_variable.py
+++ b/pybamm/expression_tree/independent_variable.py
@@ -48,7 +48,7 @@ def new_copy(self):
         """ See :meth:`pybamm.Symbol.new_copy()`. """
         return Time()
 
-    def _base_evaluate(self, t, y=None, u=None):
+    def _base_evaluate(self, t=None, y=None, y_dot=None, u=None):
         """ See :meth:`pybamm.Symbol._base_evaluate()`. """
         if t is None:
             raise ValueError("t must be provided")
diff --git a/pybamm/expression_tree/input_parameter.py b/pybamm/expression_tree/input_parameter.py
index 148c062b0c..395b84ebb2 100644
--- a/pybamm/expression_tree/input_parameter.py
+++ b/pybamm/expression_tree/input_parameter.py
@@ -36,7 +36,7 @@ def _jac(self, variable):
         """ See :meth:`pybamm.Symbol._jac()`. """
         return pybamm.Scalar(0)
 
-    def _base_evaluate(self, t=None, y=None, u=None):
+    def _base_evaluate(self, t=None, y=None, y_dot=None, u=None):
         # u should be a dictionary
         # convert 'None' to empty dictionary for more informative error
         if u is None:
diff --git a/pybamm/expression_tree/operations/convert_to_casadi.py b/pybamm/expression_tree/operations/convert_to_casadi.py
index ed8d706ee6..e76c1c330b 100644
--- a/pybamm/expression_tree/operations/convert_to_casadi.py
+++ b/pybamm/expression_tree/operations/convert_to_casadi.py
@@ -13,7 +13,7 @@ def __init__(self, casadi_symbols=None):
 
         pybamm.citations.register("Andersson2019")
 
-    def convert(self, symbol, t=None, y=None, u=None):
+    def convert(self, symbol, t, y, y_dot, u):
         """
         This function recurses down the tree, converting the PyBaMM expression tree to
         a CasADi expression tree
@@ -39,12 +39,12 @@ def convert(self, symbol, t=None, y=None, u=None):
         except KeyError:
             # Change u to empty dictionary if it's None
             u = u or {}
-            casadi_symbol = self._convert(symbol, t, y, u)
+            casadi_symbol = self._convert(symbol, t, y, y_dot, u)
             self._casadi_symbols[symbol.id] = casadi_symbol
 
             return casadi_symbol
 
-    def _convert(self, symbol, t=None, y=None, u=None):
+    def _convert(self, symbol, t, y, y_dot, u):
         """ See :meth:`CasadiConverter.convert()`. """
         if isinstance(
             symbol,
@@ -56,34 +56,41 @@ def _convert(self, symbol, t=None, y=None, u=None):
                 pybamm.ExternalVariable,
             ),
         ):
-            return casadi.MX(symbol.evaluate(t, y, u))
+            return casadi.MX(symbol.evaluate(t, y, y_dot, u))
 
         elif isinstance(symbol, pybamm.StateVector):
             if y is None:
                 raise ValueError("Must provide a 'y' for converting state vectors")
             return casadi.vertcat(*[y[y_slice] for y_slice in symbol.y_slices])
 
+        elif isinstance(symbol, pybamm.StateVectorDot):
+            if y_dot is None:
+                raise ValueError("Must provide a 'y_dot' for converting state vectors")
+            return casadi.vertcat(*[y_dot[y_slice] for y_slice in symbol.y_slices])
+
         elif isinstance(symbol, pybamm.BinaryOperator):
             left, right = symbol.children
             # process children
-            converted_left = self.convert(left, t, y, u)
-            converted_right = self.convert(right, t, y, u)
+            converted_left = self.convert(left, t, y, y_dot, u)
+            converted_right = self.convert(right, t, y, y_dot, u)
+
             if isinstance(symbol, pybamm.Minimum):
                 return casadi.fmin(converted_left, converted_right)
             if isinstance(symbol, pybamm.Maximum):
                 return casadi.fmax(converted_left, converted_right)
+
             # _binary_evaluate defined in derived classes for specific rules
             return symbol._binary_evaluate(converted_left, converted_right)
 
         elif isinstance(symbol, pybamm.UnaryOperator):
-            converted_child = self.convert(symbol.child, t, y, u)
+            converted_child = self.convert(symbol.child, t, y, y_dot, u)
             if isinstance(symbol, pybamm.AbsoluteValue):
                 return casadi.fabs(converted_child)
             return symbol._unary_evaluate(converted_child)
 
         elif isinstance(symbol, pybamm.Function):
             converted_children = [
-                self.convert(child, t, y, u) for child in symbol.children
+                self.convert(child, t, y, y_dot, u) for child in symbol.children
             ]
             # Special functions
             if symbol.function == np.min:
@@ -114,7 +121,7 @@ def _convert(self, symbol, t=None, y=None, u=None):
                 return symbol._function_evaluate(converted_children)
         elif isinstance(symbol, pybamm.Concatenation):
             converted_children = [
-                self.convert(child, t, y, u) for child in symbol.children
+                self.convert(child, t, y, y_dot, u) for child in symbol.children
             ]
             if isinstance(symbol, (pybamm.NumpyConcatenation, pybamm.SparseStack)):
                 return casadi.vertcat(*converted_children)
diff --git a/pybamm/expression_tree/operations/jacobian.py b/pybamm/expression_tree/operations/jacobian.py
index 42bd3683da..b231c0c784 100644
--- a/pybamm/expression_tree/operations/jacobian.py
+++ b/pybamm/expression_tree/operations/jacobian.py
@@ -5,8 +5,22 @@
 
 
 class Jacobian(object):
-    def __init__(self, known_jacs=None):
+    """
+    Helper class to calculate the jacobian of an expression.
+
+    Parameters
+    ----------
+
+    known_jacs: dict {variable ids -> :class:`pybamm.Symbol`}
+        cached jacobians
+
+    clear_domain: bool
+        wether or not the jacobian clears the domain (default True)
+    """
+
+    def __init__(self, known_jacs=None, clear_domain=True):
         self._known_jacs = known_jacs or {}
+        self._clear_domain = clear_domain
 
     def jac(self, symbol, variable):
         """
@@ -75,6 +89,7 @@ def _jac(self, symbol, variable):
                     )
                 )
 
-        # jacobian removes the domain(s)
-        jac.clear_domains()
+        # jacobian by default removes the domain(s)
+        if self._clear_domain:
+            jac.clear_domains()
         return jac
diff --git a/pybamm/expression_tree/scalar.py b/pybamm/expression_tree/scalar.py
index b96d618c70..3b3d1225f5 100644
--- a/pybamm/expression_tree/scalar.py
+++ b/pybamm/expression_tree/scalar.py
@@ -51,7 +51,7 @@ def set_id(self):
             (self.__class__, self.name) + tuple(self.domain) + tuple(str(self._value))
         )
 
-    def _base_evaluate(self, t=None, y=None, u=None):
+    def _base_evaluate(self, t=None, y=None, y_dot=None, u=None):
         """ See :meth:`pybamm.Symbol._base_evaluate()`. """
         return self._value
 
diff --git a/pybamm/expression_tree/state_vector.py b/pybamm/expression_tree/state_vector.py
index d99079b085..b19ee9108d 100644
--- a/pybamm/expression_tree/state_vector.py
+++ b/pybamm/expression_tree/state_vector.py
@@ -7,7 +7,7 @@
 from scipy.sparse import csr_matrix, vstack
 
 
-class StateVector(pybamm.Symbol):
+class StateVectorBase(pybamm.Symbol):
     """
     node in the expression tree that holds a slice to read from an external vector type
 
@@ -32,6 +32,7 @@ class StateVector(pybamm.Symbol):
     def __init__(
         self,
         *y_slices,
+        base_name="y",
         name=None,
         domain=None,
         auxiliary_domains=None,
@@ -42,9 +43,10 @@ def __init__(
                 raise TypeError("all y_slices must be slice objects")
         if name is None:
             if y_slices[0].start is None:
-                name = "y[:{:d}]".format(y_slice.stop)
+                name = base_name + "[:{:d}]".format(y_slice.stop)
             else:
-                name = "y[{:d}:{:d}".format(y_slices[0].start, y_slices[0].stop)
+                name = base_name + \
+                    "[{:d}:{:d}".format(y_slices[0].start, y_slices[0].stop)
             if len(y_slices) > 1:
                 name += ",{:d}:{:d}".format(y_slices[1].start, y_slices[1].stop)
                 if len(y_slices) > 2:
@@ -103,21 +105,29 @@ def set_id(self):
             + tuple(self.domain)
         )
 
-    def _base_evaluate(self, t=None, y=None, u=None):
-        """ See :meth:`pybamm.Symbol._base_evaluate()`. """
-        if y is None:
-            raise TypeError("StateVector cannot evaluate input 'y=None'")
-        if y.shape[0] < len(self.evaluation_array):
-            raise ValueError(
-                "y is too short, so value with slice is smaller than expected"
+    def _jac_diff_vector(self, variable):
+        """
+        Differentiate a slice of a StateVector of size m with respect to another slice
+        of a different StateVector of size n. This returns a (sparse) zero matrix of
+        size m x n
+
+        Parameters
+        ----------
+        variable : :class:`pybamm.Symbol`
+            The variable with respect to which to differentiate
+
+        """
+        if len(variable.y_slices) > 1:
+            raise NotImplementedError(
+                "Jacobian only implemented for a single-slice StateVector"
             )
-        else:
-            out = (y[: len(self._evaluation_array)])[self._evaluation_array]
-            if isinstance(out, np.ndarray) and out.ndim == 1:
-                out = out[:, np.newaxis]
-            return out
+        slices_size = self.y_slices[0].stop - self.y_slices[0].start
+        variable_size = variable.last_point - variable.first_point
 
-    def _jac(self, variable):
+        # Return zeros of correct size since no entries match
+        return pybamm.Matrix(csr_matrix((slices_size, variable_size)))
+
+    def _jac_same_vector(self, variable):
         """
         Differentiate a slice of a StateVector of size m with respect to another
         slice of a StateVector of size n. This returns a (sparse) matrix of size
@@ -180,3 +190,136 @@ def _evaluate_for_shape(self):
         See :meth:`pybamm.Symbol.evaluate_for_shape()`
         """
         return np.nan * np.ones((self.size, 1))
+
+
+class StateVector(StateVectorBase):
+    """
+    node in the expression tree that holds a slice to read from an external vector type
+
+    Parameters
+    ----------
+
+    y_slice: slice
+        the slice of an external y to read
+    name: str, optional
+        the name of the node
+    domain : iterable of str, optional
+        list of domains the parameter is valid over, defaults to empty list
+    auxiliary_domains : dict of str, optional
+        dictionary of auxiliary domains
+    evaluation_array : list, optional
+        List of boolean arrays representing slices. Default is None, in which case the
+        evaluation_array is computed from y_slices.
+
+    *Extends:* :class:`Array`
+    """
+
+    def __init__(
+        self,
+        *y_slices,
+        name=None,
+        domain=None,
+        auxiliary_domains=None,
+        evaluation_array=None,
+    ):
+        super().__init__(*y_slices,
+                         base_name="y", name=name, domain=domain,
+                         auxiliary_domains=auxiliary_domains,
+                         evaluation_array=evaluation_array)
+
+    def _base_evaluate(self, t=None, y=None, y_dot=None, u=None):
+        """ See :meth:`pybamm.Symbol._base_evaluate()`. """
+        if y is None:
+            raise TypeError("StateVector cannot evaluate input 'y=None'")
+        if y.shape[0] < len(self.evaluation_array):
+            raise ValueError(
+                "y is too short, so value with slice is smaller than expected"
+            )
+
+        out = (y[: len(self._evaluation_array)])[self._evaluation_array]
+        if isinstance(out, np.ndarray) and out.ndim == 1:
+            out = out[:, np.newaxis]
+        return out
+
+    def diff(self, variable):
+        if variable.id == self.id:
+            return pybamm.Scalar(1)
+        if variable.id == pybamm.t.id:
+            return StateVectorDot(*self._y_slices, name=self.name + "'",
+                                  domain=self.domain,
+                                  auxiliary_domains=self.auxiliary_domains,
+                                  evaluation_array=self.evaluation_array)
+        else:
+            return pybamm.Scalar(0)
+
+    def _jac(self, variable):
+        if isinstance(variable, pybamm.StateVector):
+            return self._jac_same_vector(variable)
+        elif isinstance(variable, pybamm.StateVectorDot):
+            return self._jac_diff_vector(variable)
+
+
+class StateVectorDot(StateVectorBase):
+    """
+    node in the expression tree that holds a slice to read from the ydot
+
+    Parameters
+    ----------
+
+    y_slice: slice
+        the slice of an external ydot to read
+    name: str, optional
+        the name of the node
+    domain : iterable of str, optional
+        list of domains the parameter is valid over, defaults to empty list
+    auxiliary_domains : dict of str, optional
+        dictionary of auxiliary domains
+    evaluation_array : list, optional
+        List of boolean arrays representing slices. Default is None, in which case the
+        evaluation_array is computed from y_slices.
+
+    *Extends:* :class:`Array`
+    """
+
+    def __init__(
+        self,
+        *y_slices,
+        name=None,
+        domain=None,
+        auxiliary_domains=None,
+        evaluation_array=None,
+    ):
+        super().__init__(*y_slices,
+                         base_name="y_dot", name=name, domain=domain,
+                         auxiliary_domains=auxiliary_domains,
+                         evaluation_array=evaluation_array)
+
+    def _base_evaluate(self, t=None, y=None, y_dot=None, u=None):
+        """ See :meth:`pybamm.Symbol._base_evaluate()`. """
+        if y_dot is None:
+            raise TypeError("StateVectorDot cannot evaluate input 'y_dot=None'")
+        if y_dot.shape[0] < len(self.evaluation_array):
+            raise ValueError(
+                "y_dot is too short, so value with slice is smaller than expected"
+            )
+
+        out = (y_dot[: len(self._evaluation_array)])[self._evaluation_array]
+        if isinstance(out, np.ndarray) and out.ndim == 1:
+            out = out[:, np.newaxis]
+        return out
+
+    def diff(self, variable):
+        if variable.id == self.id:
+            return pybamm.Scalar(1)
+        elif variable.id == pybamm.t.id:
+            raise pybamm.ModelError(
+                "cannot take second time derivative of a state vector"
+            )
+        else:
+            return pybamm.Scalar(0)
+
+    def _jac(self, variable):
+        if isinstance(variable, pybamm.StateVectorDot):
+            return self._jac_same_vector(variable)
+        elif isinstance(variable, pybamm.StateVector):
+            return self._jac_diff_vector(variable)
diff --git a/pybamm/expression_tree/symbol.py b/pybamm/expression_tree/symbol.py
index b61efc36e1..c59c69060f 100644
--- a/pybamm/expression_tree/symbol.py
+++ b/pybamm/expression_tree/symbol.py
@@ -485,6 +485,12 @@ def diff(self, variable):
             return pybamm.Scalar(1)
         elif any(variable.id == x.id for x in self.pre_order()):
             return self._diff(variable)
+        elif variable.id == pybamm.t.id and \
+                any(
+                    isinstance(x, (pybamm.VariableBase, pybamm.StateVectorBase))
+                    for x in self.pre_order()
+                ):
+            return self._diff(variable)
         else:
             return pybamm.Scalar(0)
 
@@ -492,12 +498,13 @@ def _diff(self, variable):
         "Default behaviour for differentiation, overriden by Binary and Unary Operators"
         raise NotImplementedError
 
-    def jac(self, variable, known_jacs=None):
+    def jac(self, variable, known_jacs=None, clear_domain=True):
         """
         Differentiate a symbol with respect to a (slice of) a State Vector.
         See :class:`pybamm.Jacobian`.
         """
-        return pybamm.Jacobian(known_jacs).jac(self, variable)
+        jac = pybamm.Jacobian(known_jacs, clear_domain=clear_domain)
+        return jac.jac(self, variable)
 
     def _jac(self, variable):
         """
@@ -506,7 +513,7 @@ def _jac(self, variable):
         """
         raise NotImplementedError
 
-    def _base_evaluate(self, t=None, y=None, u=None):
+    def _base_evaluate(self, t=None, y=None, y_dot=None, u=None):
         """evaluate expression tree
 
         will raise a ``NotImplementedError`` if this member function has not
@@ -520,7 +527,11 @@ def _base_evaluate(self, t=None, y=None, u=None):
             time at which to evaluate (default None)
 
         y : numpy.array, optional
-            array to evaluate when solving (default None)
+            array with state values to evaluate when solving (default None)
+
+        y_dot : numpy.array, optional
+            array with time derivatives of state values to evaluate when solving
+            (default None)
 
         """
         raise NotImplementedError(
@@ -530,7 +541,7 @@ def _base_evaluate(self, t=None, y=None, u=None):
             )
         )
 
-    def evaluate(self, t=None, y=None, u=None, known_evals=None):
+    def evaluate(self, t=None, y=None, y_dot=None, u=None, known_evals=None):
         """Evaluate expression tree (wrapper to allow using dict of known values).
         If the dict 'known_evals' is provided, the dict is searched for self.id; if
         self.id is in the keys, return that value; otherwise, evaluate using
@@ -541,7 +552,10 @@ def evaluate(self, t=None, y=None, u=None, known_evals=None):
         t : float or numeric type, optional
             time at which to evaluate (default None)
         y : numpy.array, optional
-            array to evaluate when solving (default None)
+            array with state values to evaluate when solving (default None)
+        y_dot : numpy.array, optional
+            array with time derivatives of state values to evaluate when solving
+            (default None)
         u : dict, optional
             dictionary of inputs to use when solving (default None)
         known_evals : dict, optional
@@ -556,10 +570,10 @@ def evaluate(self, t=None, y=None, u=None, known_evals=None):
         """
         if known_evals is not None:
             if self.id not in known_evals:
-                known_evals[self.id] = self._base_evaluate(t, y, u)
+                known_evals[self.id] = self._base_evaluate(t, y, y_dot, u)
             return known_evals[self.id], known_evals
         else:
-            return self._base_evaluate(t, y, u)
+            return self._base_evaluate(t, y, y_dot, u)
 
     def evaluate_for_shape(self):
         """Evaluate expression tree to find its shape. For symbols that cannot be
@@ -598,11 +612,13 @@ def is_constant(self):
         # do the search, return true if no relevent nodes are found
         return not any((isinstance(n, search_types)) for n in self.pre_order())
 
-    def evaluate_ignoring_errors(self):
+    def evaluate_ignoring_errors(self, t=0):
         """
         Evaluates the expression. If a node exists in the tree that cannot be evaluated
-        as a scalar or vector (e.g. Parameter, Variable, StateVector, InputParameter),
-        then None is returned. Otherwise the result of the evaluation is given
+        as a scalar or vector (e.g. Time, Parameter, Variable, StateVector), then None
+        is returned. If there is an InputParameter in the tree then a 1 is returned.
+        Otherwise the result of the evaluation is given.
+
 
         See Also
         --------
@@ -610,7 +626,7 @@ def evaluate_ignoring_errors(self):
 
         """
         try:
-            result = self.evaluate(t=0, u="shape test")
+            result = self.evaluate(t=t, u="shape test")
         except NotImplementedError:
             # return None if NotImplementedError is raised
             # (there is a e.g. Parameter, Variable, ... in the tree)
@@ -620,9 +636,15 @@ def evaluate_ignoring_errors(self):
             # (there is a e.g. StateVector in the tree)
             if error.args[0] == "StateVector cannot evaluate input 'y=None'":
                 return None
+            elif error.args[0] == "StateVectorDot cannot evaluate input 'y_dot=None'":
+                return None
             else:
                 raise error
         except ValueError as e:
+            # return None if specific ValueError is raised
+            # (there is a e.g. Time in the tree)
+            if e.args[0] == "t must be provided":
+                return None
             raise pybamm.ShapeError("Cannot find shape (original error: {})".format(e))
         return result
 
@@ -667,12 +689,12 @@ def simplify(self, simplified_symbols=None):
         """ Simplify the expression tree. See :class:`pybamm.Simplification`. """
         return pybamm.Simplification(simplified_symbols).simplify(self)
 
-    def to_casadi(self, t=None, y=None, u=None, casadi_symbols=None):
+    def to_casadi(self, t=None, y=None, ydot=None, u=None, casadi_symbols=None):
         """
         Convert the expression tree to a CasADi expression tree.
         See :class:`pybamm.CasadiConverter`.
         """
-        return pybamm.CasadiConverter(casadi_symbols).convert(self, t, y, u)
+        return pybamm.CasadiConverter(casadi_symbols).convert(self, t, y, ydot, u)
 
     def new_copy(self):
         """
@@ -702,7 +724,7 @@ def shape(self):
         # Try with some large y, to avoid having to use pre_order (slow)
         try:
             y = np.linspace(0.1, 0.9, int(1e4))
-            evaluated_self = self.evaluate(0, y, u="shape test")
+            evaluated_self = self.evaluate(0, y, y, u="shape test")
         # If that fails, fall back to calculating how big y should really be
         except ValueError:
             state_vectors_in_node = [
@@ -716,7 +738,7 @@ def shape(self):
                 )
                 # Pick a y that won't cause RuntimeWarnings
                 y = np.linspace(0.1, 0.9, min_y_size)
-            evaluated_self = self.evaluate(0, y)
+            evaluated_self = self.evaluate(0, y, y)
 
         # Return shape of evaluated object
         if isinstance(evaluated_self, numbers.Number):
diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py
index f8fbc9c98e..8ec8e0eb4d 100644
--- a/pybamm/expression_tree/unary_operators.py
+++ b/pybamm/expression_tree/unary_operators.py
@@ -63,15 +63,15 @@ def _unary_evaluate(self, child):
         """Perform unary operation on a child. """
         raise NotImplementedError
 
-    def evaluate(self, t=None, y=None, u=None, known_evals=None):
+    def evaluate(self, t=None, y=None, y_dot=None, u=None, known_evals=None):
         """ See :meth:`pybamm.Symbol.evaluate()`. """
         if known_evals is not None:
             if self.id not in known_evals:
-                child, known_evals = self.child.evaluate(t, y, u, known_evals)
+                child, known_evals = self.child.evaluate(t, y, y_dot, u, known_evals)
                 known_evals[self.id] = self._unary_evaluate(child)
             return known_evals[self.id], known_evals
         else:
-            child = self.child.evaluate(t, y, u)
+            child = self.child.evaluate(t, y, y_dot, u)
             return self._unary_evaluate(child)
 
     def _evaluate_for_shape(self):
diff --git a/pybamm/expression_tree/variable.py b/pybamm/expression_tree/variable.py
index 935305c924..6568b7c65e 100644
--- a/pybamm/expression_tree/variable.py
+++ b/pybamm/expression_tree/variable.py
@@ -6,7 +6,7 @@
 import numpy as np
 
 
-class Variable(pybamm.Symbol):
+class VariableBase(pybamm.Symbol):
     """A node in the expression tree represending a dependent variable
 
     This node will be discretised by :class:`.Discretisation` and converted
@@ -48,8 +48,96 @@ def _evaluate_for_shape(self):
         )
 
 
+class Variable(VariableBase):
+    """A node in the expression tree represending a dependent variable
+
+    This node will be discretised by :class:`.Discretisation` and converted
+    to a :class:`pybamm.StateVector` node.
+
+    Parameters
+    ----------
+
+    name : str
+        name of the node
+    domain : iterable of str
+        list of domains that this variable is valid over
+    auxiliary_domains : dict
+        dictionary of auxiliary domains ({'secondary': ..., 'tertiary': ...}). For
+        example, for the single particle model, the particle concentration would be a
+        Variable with domain 'negative particle' and secondary auxiliary domain 'current
+        collector'. For the DFN, the particle concentration would be a Variable with
+        domain 'negative particle', secondary domain 'negative electrode' and tertiary
+        domain 'current collector'
+
+    *Extends:* :class:`Symbol`
+    """
+
+    def __init__(self, name, domain=None, auxiliary_domains=None):
+        super().__init__(name, domain=domain, auxiliary_domains=auxiliary_domains)
+
+    def diff(self, variable):
+        if variable.id == self.id:
+            return pybamm.Scalar(1)
+        elif variable.id == pybamm.t.id:
+            return pybamm.VariableDot(self.name + "'",
+                                      domain=self.domain,
+                                      auxiliary_domains=self.auxiliary_domains)
+        else:
+            return pybamm.Scalar(0)
+
+
+class VariableDot(VariableBase):
+    """
+    A node in the expression tree represending the time derviative of a dependent
+    variable
+
+    This node will be discretised by :class:`.Discretisation` and converted
+    to a :class:`pybamm.StateVectorDot` node.
+
+    Parameters
+    ----------
+
+    name : str
+        name of the node
+    domain : iterable of str
+        list of domains that this variable is valid over
+    auxiliary_domains : dict
+        dictionary of auxiliary domains ({'secondary': ..., 'tertiary': ...}). For
+        example, for the single particle model, the particle concentration would be a
+        Variable with domain 'negative particle' and secondary auxiliary domain 'current
+        collector'. For the DFN, the particle concentration would be a Variable with
+        domain 'negative particle', secondary domain 'negative electrode' and tertiary
+        domain 'current collector'
+
+    *Extends:* :class:`Symbol`
+    """
+
+    def __init__(self, name, domain=None, auxiliary_domains=None):
+        super().__init__(name, domain=domain, auxiliary_domains=auxiliary_domains)
+
+    def get_variable(self):
+        """
+        return a :class:`.Variable` corresponding to this VariableDot
+
+        Note: Variable._jac adds a dash to the name of the corresponding VariableDot, so
+        we remove this here
+
+        """
+        return Variable(self.name[:-1],
+                        domain=self._domain,
+                        auxiliary_domains=self._auxiliary_domains)
+
+    def diff(self, variable):
+        if variable.id == self.id:
+            return pybamm.Scalar(1)
+        elif variable.id == pybamm.t.id:
+            raise pybamm.ModelError("cannot take second time derivative of a Variable")
+        else:
+            return pybamm.Scalar(0)
+
+
 class ExternalVariable(Variable):
-    """A node in the expression tree represending an external variable variable
+    """A node in the expression tree representing an external variable variable
 
     This node will be discretised by :class:`.Discretisation` and converted
     to a :class:`.Vector` node.
@@ -84,7 +172,7 @@ def _evaluate_for_shape(self):
         """ See :meth:`pybamm.Symbol.evaluate_for_shape_using_domain()` """
         return np.nan * np.ones((self.size, 1))
 
-    def _base_evaluate(self, t=None, y=None, u=None):
+    def _base_evaluate(self, t=None, y=None, y_dot=None, u=None):
         # u should be a dictionary
         # convert 'None' to empty dictionary for more informative error
         if u is None:
@@ -109,3 +197,12 @@ def _base_evaluate(self, t=None, y=None, u=None):
         # raise more informative error if can't find name in dict
         except KeyError:
             raise KeyError("External variable '{}' not found".format(self.name))
+
+    def diff(self, variable):
+        if variable.id == self.id:
+            return pybamm.Scalar(1)
+        elif variable.id == pybamm.t.id:
+            raise pybamm.ModelError(
+                "cannot take time derivative of an external variable")
+        else:
+            return pybamm.Scalar(0)
diff --git a/pybamm/models/base_model.py b/pybamm/models/base_model.py
index f8ae8c5c17..641ef347fa 100644
--- a/pybamm/models/base_model.py
+++ b/pybamm/models/base_model.py
@@ -381,6 +381,7 @@ def check_well_posedness(self, post_discretisation=False):
         post_discretisation : boolean
             A flag indicating tests to be skipped after discretisation
         """
+        self.check_for_time_derivatives()
         self.check_well_determined(post_discretisation)
         self.check_algebraic_equations(post_discretisation)
         self.check_ics_bcs()
@@ -391,6 +392,35 @@ def check_well_posedness(self, post_discretisation=False):
         if pybamm.settings.debug_mode is True and post_discretisation is False:
             self.check_variables()
 
+    def check_for_time_derivatives(self):
+        # Check that no variable time derivatives exist in the rhs equations
+        for key, eq in self.rhs.items():
+            for node in eq.pre_order():
+                if isinstance(node, pybamm.VariableDot):
+                    raise pybamm.ModelError(
+                        "time derivative of variable found ({}) in rhs equation {}"
+                        .format(node, key)
+                    )
+                if isinstance(node, pybamm.StateVectorDot):
+                    raise pybamm.ModelError(
+                        "time derivative of state vector found ({}) in rhs equation {}"
+                        .format(node, key)
+                    )
+
+        # Check that no variable time derivatives exist in the algebraic equations
+        for key, eq in self.algebraic.items():
+            for node in eq.pre_order():
+                if isinstance(node, pybamm.VariableDot):
+                    raise pybamm.ModelError(
+                        "time derivative of variable found ({}) in algebraic"
+                        "equation {}".format(node, key)
+                    )
+                if isinstance(node, pybamm.StateVectorDot):
+                    raise pybamm.ModelError(
+                        "time derivative of state vector found ({}) in algebraic"
+                        "equation {}".format(node, key)
+                    )
+
     def check_well_determined(self, post_discretisation):
         """ Check that the model is not under- or over-determined. """
         # Equations (differential and algebraic)
@@ -410,6 +440,10 @@ def check_well_determined(self, post_discretisation):
             vars_in_eqns.update(
                 [x.id for x in eqn.pre_order() if isinstance(x, pybamm.Variable)]
             )
+            vars_in_eqns.update(
+                [x.get_variable().id for x in eqn.pre_order()
+                 if isinstance(x, pybamm.VariableDot)]
+            )
         for var, eqn in self.algebraic.items():
             vars_in_algebraic_keys.update(
                 [x.id for x in var.pre_order() if isinstance(x, pybamm.Variable)]
@@ -417,11 +451,19 @@ def check_well_determined(self, post_discretisation):
             vars_in_eqns.update(
                 [x.id for x in eqn.pre_order() if isinstance(x, pybamm.Variable)]
             )
+            vars_in_eqns.update(
+                [x.get_variable().id for x in eqn.pre_order()
+                 if isinstance(x, pybamm.VariableDot)]
+            )
         for var, side_eqn in self.boundary_conditions.items():
             for side, (eqn, typ) in side_eqn.items():
                 vars_in_eqns.update(
                     [x.id for x in eqn.pre_order() if isinstance(x, pybamm.Variable)]
                 )
+                vars_in_eqns.update(
+                    [x.get_variable().id for x in eqn.pre_order()
+                     if isinstance(x, pybamm.VariableDot)]
+                )
         # If any keys are repeated between rhs and algebraic then the model is
         # overdetermined
         if not set(vars_in_rhs_keys).isdisjoint(vars_in_algebraic_keys):
diff --git a/pybamm/processed_variable.py b/pybamm/processed_variable.py
index 23d74a461c..19a3df714a 100644
--- a/pybamm/processed_variable.py
+++ b/pybamm/processed_variable.py
@@ -41,14 +41,14 @@ def __init__(self, base_variable, solution, known_evals=None):
             self.base_eval, self.known_evals[solution.t[0]] = base_variable.evaluate(
                 solution.t[0],
                 solution.y[:, 0],
-                {name: inp[0] for name, inp in solution.inputs.items()},
+                u={name: inp[0] for name, inp in solution.inputs.items()},
                 known_evals=self.known_evals[solution.t[0]],
             )
         else:
             self.base_eval = base_variable.evaluate(
                 solution.t[0],
                 solution.y[:, 0],
-                {name: inp[0] for name, inp in solution.inputs.items()},
+                u={name: inp[0] for name, inp in solution.inputs.items()},
             )
 
         # handle 2D (in space) finite element variables differently
@@ -96,10 +96,10 @@ def initialise_1D(self):
             inputs = {name: inp[idx] for name, inp in self.inputs.items()}
             if self.known_evals:
                 entries[idx], self.known_evals[t] = self.base_variable.evaluate(
-                    t, u, inputs, known_evals=self.known_evals[t]
+                    t, u, u=inputs, known_evals=self.known_evals[t]
                 )
             else:
-                entries[idx] = self.base_variable.evaluate(t, u, inputs)
+                entries[idx] = self.base_variable.evaluate(t, u, u=inputs)
 
         # No discretisation provided, or variable has no domain (function of t only)
         self._interpolation_function = interp.interp1d(
@@ -120,12 +120,12 @@ def initialise_2D(self):
             inputs = {name: inp[idx] for name, inp in self.inputs.items()}
             if self.known_evals:
                 eval_and_known_evals = self.base_variable.evaluate(
-                    t, u, inputs, known_evals=self.known_evals[t]
+                    t, u, u=inputs, known_evals=self.known_evals[t]
                 )
                 entries[:, idx] = eval_and_known_evals[0][:, 0]
                 self.known_evals[t] = eval_and_known_evals[1]
             else:
-                entries[:, idx] = self.base_variable.evaluate(t, u, inputs)[:, 0]
+                entries[:, idx] = self.base_variable.evaluate(t, u, u=inputs)[:, 0]
 
         # Process the discretisation to get x values
         nodes = self.mesh[0].nodes
@@ -226,7 +226,7 @@ def initialise_3D(self):
             inputs = {name: inp[idx] for name, inp in self.inputs.items()}
             if self.known_evals:
                 eval_and_known_evals = self.base_variable.evaluate(
-                    t, u, inputs, known_evals=self.known_evals[t]
+                    t, u, u=inputs, known_evals=self.known_evals[t]
                 )
                 entries[:, :, idx] = np.reshape(
                     eval_and_known_evals[0],
@@ -236,7 +236,7 @@ def initialise_3D(self):
                 self.known_evals[t] = eval_and_known_evals[1]
             else:
                 entries[:, :, idx] = np.reshape(
-                    self.base_variable.evaluate(t, u, inputs),
+                    self.base_variable.evaluate(t, u, u=inputs),
                     [first_dim_size, second_dim_size],
                     order="F",
                 )
@@ -265,7 +265,7 @@ def initialise_2Dspace_scikit_fem(self):
         inputs = {name: inp[0] for name, inp in self.inputs.items()}
 
         entries = np.reshape(
-            self.base_variable.evaluate(0, self.u_sol, inputs), [len_y, len_z]
+            self.base_variable.evaluate(0, self.u_sol, u=inputs), [len_y, len_z]
         )
 
         # assign attributes for reference
@@ -296,13 +296,13 @@ def initialise_3D_scikit_fem(self):
 
             if self.known_evals:
                 eval_and_known_evals = self.base_variable.evaluate(
-                    t, u, inputs, known_evals=self.known_evals[t]
+                    t, u, u=inputs, known_evals=self.known_evals[t]
                 )
                 entries[:, :, idx] = np.reshape(eval_and_known_evals[0], [len_y, len_z])
                 self.known_evals[t] = eval_and_known_evals[1]
             else:
                 entries[:, :, idx] = np.reshape(
-                    self.base_variable.evaluate(t, u, inputs), [len_y, len_z]
+                    self.base_variable.evaluate(t, u, u=inputs), [len_y, len_z]
                 )
 
         # assign attributes for reference
diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index f849d8b3c5..f797909580 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -130,7 +130,7 @@ def set_up(self, model, inputs=None):
             )
 
         inputs = inputs or {}
-        y0 = model.concatenated_initial_conditions.evaluate(0, None, inputs)
+        y0 = model.concatenated_initial_conditions.evaluate(0, None, u=inputs)
 
         # Set model timescale
         model.timescale_eval = model.timescale.evaluate(u=inputs)
@@ -164,10 +164,10 @@ def set_up(self, model, inputs=None):
             # Convert model attributes to casadi
             t_casadi = casadi.MX.sym("t")
             y_diff = casadi.MX.sym(
-                "y_diff", len(model.concatenated_rhs.evaluate(0, y0, inputs))
+                "y_diff", len(model.concatenated_rhs.evaluate(0, y0, u=inputs))
             )
             y_alg = casadi.MX.sym(
-                "y_alg", len(model.concatenated_algebraic.evaluate(0, y0, inputs))
+                "y_alg", len(model.concatenated_algebraic.evaluate(0, y0, u=inputs))
             )
             y_casadi = casadi.vertcat(y_diff, y_alg)
             u_casadi = {}
@@ -210,7 +210,7 @@ def report(string):
             else:
                 # Process with CasADi
                 report(f"Converting {name} to CasADi")
-                func = func.to_casadi(t_casadi, y_casadi, u_casadi)
+                func = func.to_casadi(t_casadi, y_casadi, u=u_casadi)
                 if use_jacobian:
                     report(f"Calculating jacobian for {name} using CasADi")
                     jac_casadi = casadi.jacobian(func, y_casadi)
@@ -810,7 +810,7 @@ def get_termination_reason(self, solution, events):
                         event.expression.evaluate(
                             solution.t_event,
                             solution.y_event,
-                            {k: v[-1] for k, v in solution.inputs.items()},
+                            u={k: v[-1] for k, v in solution.inputs.items()},
                         )
                     )
             termination_event = min(final_event_values, key=final_event_values.get)
@@ -855,11 +855,12 @@ def __call__(self, t, y):
 
     def function(self, t, y):
         if self.form == "casadi":
+            states_eval = self._function(t, y, self.inputs)
             if self.name in ["RHS", "algebraic", "residuals", "event"]:
-                return self._function(t, y, self.inputs).full()
+                return states_eval.full()
             else:
                 # keep jacobians sparse
-                return self._function(t, y, self.inputs)
+                return states_eval
         else:
             return self._function(t, y, self.inputs, known_evals={})[0]
 
diff --git a/tests/integration/test_solvers/test_external_variables.py b/tests/integration/test_solvers/test_external_variables.py
index cd2a43abb9..f8f796284f 100644
--- a/tests/integration/test_solvers/test_external_variables.py
+++ b/tests/integration/test_solvers/test_external_variables.py
@@ -50,7 +50,7 @@ def test_external_variables_SPMe(self):
         t = sim.solution.t[-1]
         y = sim.solution.y[:, -1]
         u = external_variables
-        sim.built_model.variables[var].evaluate(t, y, u)
+        sim.built_model.variables[var].evaluate(t, y, u=u)
         sim.solution[var](t)
 
 
diff --git a/tests/unit/test_discretisations/test_discretisation.py b/tests/unit/test_discretisations/test_discretisation.py
index da38f07782..cb0f015a99 100644
--- a/tests/unit/test_discretisations/test_discretisation.py
+++ b/tests/unit/test_discretisations/test_discretisation.py
@@ -329,6 +329,13 @@ def test_process_symbol_base(self):
         var_disc = disc.process_symbol(var)
         self.assertIsInstance(var_disc, pybamm.StateVector)
         self.assertEqual(var_disc.y_slices[0], disc.y_slices[var.id][0])
+
+        # variable dot
+        var_dot = pybamm.VariableDot("var'")
+        var_dot_disc = disc.process_symbol(var_dot)
+        self.assertIsInstance(var_dot_disc, pybamm.StateVectorDot)
+        self.assertEqual(var_dot_disc.y_slices[0], disc.y_slices[var.id][0])
+
         # scalar
         scal = pybamm.Scalar(5)
         scal_disc = disc.process_symbol(scal)
@@ -693,6 +700,25 @@ def test_process_model_ode(self):
         with self.assertRaises(pybamm.ModelError):
             disc.process_model(model)
 
+        # test that any time derivatives of variables in rhs raises an
+        # error
+        model = pybamm.BaseModel()
+        model.rhs = {c: pybamm.div(N) + c.diff(pybamm.t),
+                     T: pybamm.div(q), S: pybamm.div(p)}
+        model.initial_conditions = {
+            c: pybamm.Scalar(2),
+            T: pybamm.Scalar(5),
+            S: pybamm.Scalar(8),
+        }
+        model.boundary_conditions = {
+            c: {"left": (0, "Neumann"), "right": (0, "Neumann")},
+            T: {"left": (0, "Neumann"), "right": (0, "Neumann")},
+            S: {"left": (0, "Neumann"), "right": (0, "Neumann")},
+        }
+        model.variables = {"ST": S * T}
+        with self.assertRaises(pybamm.ModelError):
+            disc.process_model(model)
+
     def test_process_model_fail(self):
         # one equation
         c = pybamm.Variable("c")
@@ -807,6 +833,20 @@ def test_process_model_dae(self):
         jacobian = expr.evaluate(0, y0, known_evals=known_evals)[0]
         np.testing.assert_array_equal(jacobian_actual, jacobian.toarray())
 
+        # check that any time derivatives of variables in algebraic raises an
+        # error
+        model = pybamm.BaseModel()
+        model.rhs = {c: pybamm.div(N)}
+        model.algebraic = {d: d - 2 * c.diff(pybamm.t)}
+        model.initial_conditions = {d: pybamm.Scalar(6), c: pybamm.Scalar(3)}
+        model.boundary_conditions = {
+            c: {"left": (0, "Neumann"), "right": (0, "Neumann")}
+        }
+        model.variables = {"c": c, "N": N, "d": d}
+
+        with self.assertRaises(pybamm.ModelError):
+            disc.process_model(model)
+
     def test_process_model_concatenation(self):
         # concatenation of variables as the key
         cn = pybamm.Variable("c", domain=["negative electrode"])
diff --git a/tests/unit/test_expression_tree/test_d_dt.py b/tests/unit/test_expression_tree/test_d_dt.py
new file mode 100644
index 0000000000..f9eb5a2cd9
--- /dev/null
+++ b/tests/unit/test_expression_tree/test_d_dt.py
@@ -0,0 +1,70 @@
+#
+# Tests for the Scalar class
+#
+import pybamm
+
+import unittest
+import numpy as np
+
+
+class TestDDT(unittest.TestCase):
+    def test_time_derivative(self):
+        a = pybamm.Scalar(5).diff(pybamm.t)
+        self.assertIsInstance(a, pybamm.Scalar)
+        self.assertEqual(a.value, 0)
+
+        a = pybamm.t.diff(pybamm.t)
+        self.assertIsInstance(a, pybamm.Scalar)
+        self.assertEqual(a.value, 1)
+
+        a = (pybamm.t**2).diff(pybamm.t)
+        self.assertEqual(a.id, (2 * pybamm.t ** 1 * 1).id)
+        self.assertEqual(a.simplify().id, (2 * pybamm.t).id)
+        self.assertEqual(a.evaluate(t=1), 2)
+
+        a = (2 + pybamm.t**2).diff(pybamm.t)
+        self.assertEqual(a.simplify().id, (2 * pybamm.t).id)
+        self.assertEqual(a.evaluate(t=1), 2)
+
+    def test_time_derivative_of_variable(self):
+
+        a = (pybamm.Variable('a')).diff(pybamm.t)
+        self.assertIsInstance(a, pybamm.VariableDot)
+        self.assertEqual(a.name, "a'")
+
+        p = pybamm.Parameter('p')
+        a = (1 + p * pybamm.Variable('a')).diff(pybamm.t).simplify()
+        self.assertIsInstance(a, pybamm.Multiplication)
+        self.assertEqual(a.children[0].name, 'p')
+        self.assertEqual(a.children[1].name, "a'")
+
+        with self.assertRaises(pybamm.ModelError):
+            a = (pybamm.Variable('a')).diff(pybamm.t).diff(pybamm.t)
+
+        with self.assertRaises(pybamm.ModelError):
+            a = pybamm.ExternalVariable("a", 1).diff(pybamm.t)
+
+    def test_time_derivative_of_state_vector(self):
+
+        sv = pybamm.StateVector(slice(0, 10))
+        y_dot = np.linspace(0, 2, 19)
+
+        a = sv.diff(pybamm.t)
+        self.assertIsInstance(a, pybamm.StateVectorDot)
+        self.assertEqual(a.name[-1], "'")
+        np.testing.assert_array_equal(
+            a.evaluate(y_dot=y_dot), np.linspace(0, 1, 10)[:, np.newaxis]
+        )
+
+        with self.assertRaises(pybamm.ModelError):
+            a = (sv).diff(pybamm.t).diff(pybamm.t)
+
+
+if __name__ == "__main__":
+    print("Add -v for more debug output")
+    import sys
+
+    if "-v" in sys.argv:
+        debug = True
+    pybamm.settings.debug_mode = True
+    unittest.main()
diff --git a/tests/unit/test_expression_tree/test_operations/test_convert_to_casadi.py b/tests/unit/test_expression_tree/test_operations/test_convert_to_casadi.py
index 915a449c31..fa9992ad24 100644
--- a/tests/unit/test_expression_tree/test_operations/test_convert_to_casadi.py
+++ b/tests/unit/test_expression_tree/test_operations/test_convert_to_casadi.py
@@ -66,9 +66,11 @@ def test_convert_array_symbols(self):
 
         casadi_t = casadi.MX.sym("t")
         casadi_y = casadi.MX.sym("y", 10)
+        casadi_y_dot = casadi.MX.sym("y_dot", 10)
 
         pybamm_t = pybamm.Time()
         pybamm_y = pybamm.StateVector(slice(0, 10))
+        pybamm_y_dot = pybamm.StateVectorDot(slice(0, 10))
 
         # Time
         self.assertEqual(pybamm_t.to_casadi(casadi_t, casadi_y), casadi_t)
@@ -76,6 +78,12 @@ def test_convert_array_symbols(self):
         # State Vector
         self.assert_casadi_equal(pybamm_y.to_casadi(casadi_t, casadi_y), casadi_y)
 
+        # State Vector Dot
+        self.assert_casadi_equal(
+            pybamm_y_dot.to_casadi(casadi_t, casadi_y, casadi_y_dot),
+            casadi_y_dot
+        )
+
     def test_special_functions(self):
         a = pybamm.Array(np.array([1, 2, 3, 4, 5]))
         self.assert_casadi_equal(pybamm.max(a).to_casadi(), casadi.MX(5), evalf=True)
@@ -159,6 +167,7 @@ def myfunction(x, y):
     def test_convert_input_parameter(self):
         casadi_t = casadi.MX.sym("t")
         casadi_y = casadi.MX.sym("y", 10)
+        casadi_ydot = casadi.MX.sym("ydot", 10)
         casadi_us = {
             "Input 1": casadi.MX.sym("Input 1"),
             "Input 2": casadi.MX.sym("Input 2"),
@@ -170,18 +179,19 @@ def test_convert_input_parameter(self):
 
         # Input only
         self.assert_casadi_equal(
-            pybamm_u1.to_casadi(casadi_t, casadi_y, casadi_us), casadi_us["Input 1"]
+            pybamm_u1.to_casadi(casadi_t, casadi_y, casadi_ydot, casadi_us),
+            casadi_us["Input 1"]
         )
 
         # More complex
         expr = pybamm_u1 + pybamm_y
         self.assert_casadi_equal(
-            expr.to_casadi(casadi_t, casadi_y, casadi_us),
+            expr.to_casadi(casadi_t, casadi_y, casadi_ydot, casadi_us),
             casadi_us["Input 1"] + casadi_y,
         )
         expr = pybamm_u2 * pybamm_y
         self.assert_casadi_equal(
-            expr.to_casadi(casadi_t, casadi_y, casadi_us),
+            expr.to_casadi(casadi_t, casadi_y, casadi_ydot, casadi_us),
             casadi_us["Input 2"] * casadi_y,
         )
 
@@ -199,13 +209,14 @@ def test_convert_external_variable(self):
 
         # External only
         self.assert_casadi_equal(
-            pybamm_u1.to_casadi(casadi_t, casadi_y, casadi_us), casadi_us["External 1"]
+            pybamm_u1.to_casadi(casadi_t, casadi_y, u=casadi_us),
+            casadi_us["External 1"]
         )
 
         # More complex
         expr = pybamm_u2 + pybamm_y
         self.assert_casadi_equal(
-            expr.to_casadi(casadi_t, casadi_y, casadi_us),
+            expr.to_casadi(casadi_t, casadi_y, u=casadi_us),
             casadi_us["External 2"] + casadi_y,
         )
 
@@ -215,6 +226,11 @@ def test_errors(self):
             ValueError, "Must provide a 'y' for converting state vectors"
         ):
             y.to_casadi()
+        y_dot = pybamm.StateVectorDot(slice(0, 10))
+        with self.assertRaisesRegex(
+            ValueError, "Must provide a 'y_dot' for converting state vectors"
+        ):
+            y_dot.to_casadi()
         var = pybamm.Variable("var")
         with self.assertRaisesRegex(TypeError, "Cannot convert symbol of type"):
             var.to_casadi()
diff --git a/tests/unit/test_expression_tree/test_operations/test_jac.py b/tests/unit/test_expression_tree/test_operations/test_jac.py
index 2711e06d29..56895d26b4 100644
--- a/tests/unit/test_expression_tree/test_operations/test_jac.py
+++ b/tests/unit/test_expression_tree/test_operations/test_jac.py
@@ -105,6 +105,58 @@ def test_nonlinear(self):
         dfunc_dy = func.jac(y).evaluate(y=y0)
         np.testing.assert_array_equal(jacobian, dfunc_dy.toarray())
 
+    def test_multislice_raises(self):
+        y1 = pybamm.StateVector(slice(0, 4), slice(7, 8))
+        y_dot1 = pybamm.StateVectorDot(slice(0, 4), slice(7, 8))
+        y2 = pybamm.StateVector(slice(4, 7))
+        with self.assertRaises(NotImplementedError):
+            y1.jac(y1)
+        with self.assertRaises(NotImplementedError):
+            y2.jac(y1)
+        with self.assertRaises(NotImplementedError):
+            y_dot1.jac(y1)
+
+    def test_linear_ydot(self):
+        y = pybamm.StateVector(slice(0, 4))
+        y_dot = pybamm.StateVectorDot(slice(0, 4))
+        u = pybamm.StateVector(slice(0, 2))
+        v = pybamm.StateVector(slice(2, 4))
+        u_dot = pybamm.StateVectorDot(slice(0, 2))
+        v_dot = pybamm.StateVectorDot(slice(2, 4))
+
+        y0 = np.ones(4)
+        y_dot0 = np.ones(4)
+
+        func = u_dot
+        jacobian = np.array([[1, 0, 0, 0], [0, 1, 0, 0]])
+        dfunc_dy = func.jac(y_dot).evaluate(y=y0, y_dot=y_dot0)
+        np.testing.assert_array_equal(jacobian, dfunc_dy.toarray())
+
+        func = -v_dot
+        jacobian = np.array([[0, 0, -1, 0], [0, 0, 0, -1]])
+        dfunc_dy = func.jac(y_dot).evaluate(y=y0, y_dot=y_dot0)
+        np.testing.assert_array_equal(jacobian, dfunc_dy.toarray())
+
+        func = u_dot
+        jacobian = np.array([[0, 0, 0, 0], [0, 0, 0, 0]])
+        dfunc_dy = func.jac(y).evaluate(y=y0, y_dot=y_dot0)
+        np.testing.assert_array_equal(jacobian, dfunc_dy.toarray())
+
+        func = -v_dot
+        jacobian = np.array([[0, 0, 0, 0], [0, 0, 0, 0]])
+        dfunc_dy = func.jac(y).evaluate(y=y0, y_dot=y_dot0)
+        np.testing.assert_array_equal(jacobian, dfunc_dy.toarray())
+
+        func = u
+        jacobian = np.array([[0, 0, 0, 0], [0, 0, 0, 0]])
+        dfunc_dy = func.jac(y_dot).evaluate(y=y0, y_dot=y_dot0)
+        np.testing.assert_array_equal(jacobian, dfunc_dy.toarray())
+
+        func = -v
+        jacobian = np.array([[0, 0, 0, 0], [0, 0, 0, 0]])
+        dfunc_dy = func.jac(y_dot).evaluate(y=y0, y_dot=y_dot0)
+        np.testing.assert_array_equal(jacobian, dfunc_dy.toarray())
+
     def test_functions(self):
         y = pybamm.StateVector(slice(0, 4))
         u = pybamm.StateVector(slice(0, 2))
diff --git a/tests/unit/test_expression_tree/test_operations/test_simplify.py b/tests/unit/test_expression_tree/test_operations/test_simplify.py
index 129aef5e98..2cd37f7537 100644
--- a/tests/unit/test_expression_tree/test_operations/test_simplify.py
+++ b/tests/unit/test_expression_tree/test_operations/test_simplify.py
@@ -16,6 +16,7 @@ def test_symbol_simplify(self):
         d = pybamm.Scalar(-1)
         e = pybamm.Scalar(2)
         g = pybamm.Variable("g")
+        gdot = pybamm.VariableDot("g'")
 
         # negate
         self.assertIsInstance((-a).simplify(), pybamm.Scalar)
@@ -175,6 +176,18 @@ def myfunction(x, y):
         self.assertIsInstance(expr.children[1], pybamm.Negate)
         self.assertIsInstance(expr.children[1].children[0], pybamm.Parameter)
 
+        expr = (e * g * b).simplify()
+        self.assertIsInstance(expr, pybamm.Multiplication)
+        self.assertIsInstance(expr.children[0], pybamm.Scalar)
+        self.assertEqual(expr.children[0].evaluate(), 2.0)
+        self.assertIsInstance(expr.children[1], pybamm.Variable)
+
+        expr = (e * gdot * b).simplify()
+        self.assertIsInstance(expr, pybamm.Multiplication)
+        self.assertIsInstance(expr.children[0], pybamm.Scalar)
+        self.assertEqual(expr.children[0].evaluate(), 2.0)
+        self.assertIsInstance(expr.children[1], pybamm.VariableDot)
+
         expr = (e + (g - c)).simplify()
         self.assertIsInstance(expr, pybamm.Addition)
         self.assertIsInstance(expr.children[0], pybamm.Scalar)
diff --git a/tests/unit/test_expression_tree/test_state_vector.py b/tests/unit/test_expression_tree/test_state_vector.py
index f6809565d3..38f3cd2c39 100644
--- a/tests/unit/test_expression_tree/test_state_vector.py
+++ b/tests/unit/test_expression_tree/test_state_vector.py
@@ -62,6 +62,27 @@ def test_failure(self):
             pybamm.StateVector(slice(0, 10), 1)
 
 
+class TestStateVectorDot(unittest.TestCase):
+    def test_evaluate(self):
+        sv = pybamm.StateVectorDot(slice(0, 10))
+        y_dot = np.linspace(0, 2, 19)
+        np.testing.assert_array_equal(
+            sv.evaluate(y_dot=y_dot), np.linspace(0, 1, 10)[:, np.newaxis]
+        )
+
+        # Try evaluating with a y that is too short
+        y_dot2 = np.ones(5)
+        with self.assertRaisesRegex(
+            ValueError,
+            "y_dot is too short, so value with slice is smaller than expected"
+        ):
+            sv.evaluate(y_dot=y_dot2)
+
+    def test_name(self):
+        sv = pybamm.StateVectorDot(slice(0, 10))
+        self.assertEqual(sv.name, "y_dot[0:10]")
+
+
 if __name__ == "__main__":
     print("Add -v for more debug output")
     import sys
diff --git a/tests/unit/test_expression_tree/test_symbol.py b/tests/unit/test_expression_tree/test_symbol.py
index dbff68dd3a..7faa4992c4 100644
--- a/tests/unit/test_expression_tree/test_symbol.py
+++ b/tests/unit/test_expression_tree/test_symbol.py
@@ -168,6 +168,13 @@ def test_symbol_evaluation(self):
         with self.assertRaises(NotImplementedError):
             a.evaluate()
 
+    def test_evaluate_ignoring_errors(self):
+        self.assertIsNone(pybamm.t.evaluate_ignoring_errors(t=None))
+        self.assertEqual(pybamm.t.evaluate_ignoring_errors(t=0), 0)
+        self.assertIsNone(pybamm.Parameter("a").evaluate_ignoring_errors())
+        self.assertIsNone(pybamm.StateVector(slice(0, 1)).evaluate_ignoring_errors())
+        self.assertEqual(pybamm.InputParameter("a").evaluate_ignoring_errors(), 1)
+
     def test_symbol_is_constant(self):
         a = pybamm.Variable("a")
         self.assertFalse(a.is_constant())
diff --git a/tests/unit/test_expression_tree/test_symbolic_diff.py b/tests/unit/test_expression_tree/test_symbolic_diff.py
index 68370ddfd3..25c6b76749 100644
--- a/tests/unit/test_expression_tree/test_symbolic_diff.py
+++ b/tests/unit/test_expression_tree/test_symbolic_diff.py
@@ -60,6 +60,12 @@ def test_diff_zero(self):
         self.assertEqual(func.diff(b).id, pybamm.Scalar(0).id)
         self.assertNotEqual(func.diff(a).id, pybamm.Scalar(0).id)
 
+    def test_diff_state_vector_dot(self):
+        a = pybamm.StateVectorDot(slice(0, 1))
+        b = pybamm.StateVector(slice(1, 2))
+        self.assertEqual(a.diff(a).id, pybamm.Scalar(1).id)
+        self.assertEqual(a.diff(b).id, pybamm.Scalar(0).id)
+
     def test_diff_heaviside(self):
         a = pybamm.Scalar(1)
         b = pybamm.StateVector(slice(0, 1))
diff --git a/tests/unit/test_expression_tree/test_variable.py b/tests/unit/test_expression_tree/test_variable.py
index 1e4e1b0bee..33f83f02b2 100644
--- a/tests/unit/test_expression_tree/test_variable.py
+++ b/tests/unit/test_expression_tree/test_variable.py
@@ -16,6 +16,14 @@ def test_variable_init(self):
         self.assertEqual(a.domain[0], "test")
         self.assertRaises(TypeError, pybamm.Variable("a", domain="test"))
 
+    def test_variable_diff(self):
+        a = pybamm.Variable("a")
+        b = pybamm.Variable("b")
+        self.assertIsInstance(a.diff(a), pybamm.Scalar)
+        self.assertEqual(a.diff(a).evaluate(), 1)
+        self.assertIsInstance(a.diff(b), pybamm.Scalar)
+        self.assertEqual(a.diff(b).evaluate(), 0)
+
     def test_variable_id(self):
         a1 = pybamm.Variable("a", domain=["negative electrode"])
         a2 = pybamm.Variable("a", domain=["negative electrode"])
@@ -26,6 +34,33 @@ def test_variable_id(self):
         self.assertNotEqual(a1.id, a4.id)
 
 
+class TestVariableDot(unittest.TestCase):
+    def test_variable_init(self):
+        a = pybamm.VariableDot("a'")
+        self.assertEqual(a.name, "a'")
+        self.assertEqual(a.domain, [])
+        a = pybamm.VariableDot("a", domain=["test"])
+        self.assertEqual(a.domain[0], "test")
+        self.assertRaises(TypeError, pybamm.Variable("a", domain="test"))
+
+    def test_variable_id(self):
+        a1 = pybamm.VariableDot("a", domain=["negative electrode"])
+        a2 = pybamm.VariableDot("a", domain=["negative electrode"])
+        self.assertEqual(a1.id, a2.id)
+        a3 = pybamm.VariableDot("b", domain=["negative electrode"])
+        a4 = pybamm.VariableDot("a", domain=["positive electrode"])
+        self.assertNotEqual(a1.id, a3.id)
+        self.assertNotEqual(a1.id, a4.id)
+
+    def test_variable_diff(self):
+        a = pybamm.VariableDot("a")
+        b = pybamm.Variable("b")
+        self.assertIsInstance(a.diff(a), pybamm.Scalar)
+        self.assertEqual(a.diff(a).evaluate(), 1)
+        self.assertIsInstance(a.diff(b), pybamm.Scalar)
+        self.assertEqual(a.diff(b).evaluate(), 0)
+
+
 class TestExternalVariable(unittest.TestCase):
     def test_external_variable_scalar(self):
         a = pybamm.ExternalVariable("a", 1)
@@ -50,6 +85,14 @@ def test_external_variable_vector(self):
         with self.assertRaisesRegex(ValueError, "External variable"):
             a.evaluate(u={"a": np.ones((5, 1))})
 
+    def test_external_variable_diff(self):
+        a = pybamm.ExternalVariable("a", 10)
+        b = pybamm.Variable("b")
+        self.assertIsInstance(a.diff(a), pybamm.Scalar)
+        self.assertEqual(a.diff(a).evaluate(), 1)
+        self.assertIsInstance(a.diff(b), pybamm.Scalar)
+        self.assertEqual(a.diff(b).evaluate(), 0)
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/unit/test_models/test_base_model.py b/tests/unit/test_models/test_base_model.py
index b9c5e6c703..2f2b37d05a 100644
--- a/tests/unit/test_models/test_base_model.py
+++ b/tests/unit/test_models/test_base_model.py
@@ -248,6 +248,61 @@ def test_check_well_posedness_variables(self):
         ):
             model.check_well_posedness(post_discretisation=True)
 
+        # model must be in semi-explicit form
+        model = pybamm.BaseModel()
+        model.rhs = {c: d.diff(pybamm.t), d: -1}
+        model.initial_conditions = {c: 1, d: 1}
+        with self.assertRaisesRegex(
+            pybamm.ModelError,
+            "time derivative of variable found",
+        ):
+            model.check_well_posedness()
+
+        # model must be in semi-explicit form
+        model = pybamm.BaseModel()
+        model.algebraic = {
+            c: 2 * d - c,
+            d: c * d.diff(pybamm.t) - d,
+        }
+        model.initial_conditions = {c: 1, d: 1}
+        with self.assertRaisesRegex(
+            pybamm.ModelError,
+            "time derivative of variable found",
+        ):
+            model.check_well_posedness()
+
+        # model must be in semi-explicit form
+        model = pybamm.BaseModel()
+        model.rhs = {c: d.diff(pybamm.t), d: -1}
+        model.initial_conditions = {c: 1, d: 1}
+        with self.assertRaisesRegex(
+            pybamm.ModelError,
+            "time derivative of variable found",
+        ):
+            model.check_well_posedness()
+
+        # model must be in semi-explicit form
+        model = pybamm.BaseModel()
+        model.algebraic = {
+            d: 5 * pybamm.StateVector(slice(0, 15)) - 1,
+            c: 5 * pybamm.StateVectorDot(slice(0, 15)) - 1
+        }
+        with self.assertRaisesRegex(
+            pybamm.ModelError,
+            "time derivative of state vector found",
+        ):
+            model.check_well_posedness(post_discretisation=True)
+
+        # model must be in semi-explicit form
+        model = pybamm.BaseModel()
+        model.rhs = {c: 5 * pybamm.StateVectorDot(slice(0, 15)) - 1}
+        model.initial_conditions = {c: 1}
+        with self.assertRaisesRegex(
+            pybamm.ModelError,
+            "time derivative of state vector found",
+        ):
+            model.check_well_posedness(post_discretisation=True)
+
     def test_check_well_posedness_initial_boundary_conditions(self):
         # Well-posed model - Dirichlet
         whole_cell = ["negative electrode", "separator", "positive electrode"]