diff --git a/CHANGELOG.md b/CHANGELOG.md
index 324461980a..d1e646d411 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,7 @@
 
 ## Features
 
+-   Added functionality to solver to automatically discretise a 0D model ([#947](https://github.com/pybamm-team/PyBaMM/pull/947))
 -   Made `QuickPlot` compatible with Google Colab ([#935](https://github.com/pybamm-team/PyBaMM/pull/935))
 -   Added `BasicFull` model for lead-acid ([#932](https://github.com/pybamm-team/PyBaMM/pull/932))
 
diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index 9a84cffd3e..3612df32b3 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -133,6 +133,13 @@ def process_model(self, model, inplace=True, check_model=True):
         # Prepare discretisation
         # set variables (we require the full variable not just id)
         variables = list(model.rhs.keys()) + list(model.algebraic.keys())
+        if self.spatial_methods == {} and any(var.domain != [] for var in variables):
+            for var in variables:
+                if var.domain != []:
+                    raise pybamm.DiscretisationError(
+                        "Spatial method has not been given "
+                        "for variable {} with domain {}".format(var.name, var.domain)
+                    )
 
         # Set the y split for variables
         pybamm.logger.info("Set variable slices for {}".format(model.name))
@@ -203,6 +210,9 @@ def process_model(self, model, inplace=True, check_model=True):
 
         pybamm.logger.info("Finish discretising {}".format(model.name))
 
+        # Record that the model has been discretised
+        model_disc.is_discretised = True
+
         return model_disc
 
     def set_variable_slices(self, variables):
diff --git a/pybamm/expression_tree/exceptions.py b/pybamm/expression_tree/exceptions.py
index d0f4341a20..41995623c6 100644
--- a/pybamm/expression_tree/exceptions.py
+++ b/pybamm/expression_tree/exceptions.py
@@ -67,3 +67,11 @@ class InputError(Exception):
     """
 
     pass
+
+
+class DiscretisationError(Exception):
+    """
+    A model could not be discretised
+    """
+
+    pass
diff --git a/pybamm/models/base_model.py b/pybamm/models/base_model.py
index 84776f956c..3e6fdd9b7d 100644
--- a/pybamm/models/base_model.py
+++ b/pybamm/models/base_model.py
@@ -121,6 +121,9 @@ def __init__(self, name="Unnamed model"):
         self.use_simplify = True
         self.convert_to_format = "casadi"
 
+        # Model is not initially discretised
+        self.is_discretised = False
+
         # Default timescale is 1 second
         self.timescale = pybamm.Scalar(1)
 
diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 7549f887ba..e32c141045 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -128,6 +128,18 @@ def set_up(self, model, inputs=None):
             raise pybamm.SolverError(
                 """Cannot use algebraic solver to solve model with time derivatives"""
             )
+        # Discretise model if it isn't already discretised
+        # This only works with purely 0D models, as otherwise the mesh and spatial
+        # method should be specified by the user
+        if model.is_discretised is False:
+            try:
+                disc = pybamm.Discretisation()
+                disc.process_model(model)
+            except pybamm.DiscretisationError as e:
+                raise pybamm.DiscretisationError(
+                    "Cannot automatically discretise model, "
+                    "model should be discretised before solving ({})".format(e)
+                )
 
         inputs = inputs or {}
         y0 = model.concatenated_initial_conditions.evaluate(0, None, inputs=inputs)
@@ -564,11 +576,7 @@ def solve(self, model, t_eval=None, external_variables=None, inputs=None):
         ]
 
         # remove any discontinuities after end of t_eval
-        discontinuities = [
-            v
-            for v in discontinuities
-            if v < t_eval_dimensionless[-1]
-        ]
+        discontinuities = [v for v in discontinuities if v < t_eval_dimensionless[-1]]
 
         if len(discontinuities) > 0:
             pybamm.logger.info(
diff --git a/tests/unit/test_solvers/test_base_solver.py b/tests/unit/test_solvers/test_base_solver.py
index ab051a9d06..3a7c2f130c 100644
--- a/tests/unit/test_solvers/test_base_solver.py
+++ b/tests/unit/test_solvers/test_base_solver.py
@@ -179,6 +179,29 @@ def algebraic_eval(self, t, y, inputs):
         ):
             solver.calculate_consistent_state(Model())
 
+    def test_discretise_model(self):
+        # Make sure 0D model is automatically discretised
+        model = pybamm.BaseModel()
+        v = pybamm.Variable("v")
+        model.rhs = {v: -1}
+        model.initial_conditions = {v: 1}
+
+        solver = pybamm.BaseSolver()
+        self.assertFalse(model.is_discretised)
+        solver.set_up(model, {})
+        self.assertTrue(model.is_discretised)
+
+        # 1D model cannot be automatically discretised
+        model = pybamm.BaseModel()
+        v = pybamm.Variable("v", domain="line")
+        model.rhs = {v: -1}
+        model.initial_conditions = {v: 1}
+
+        with self.assertRaisesRegex(
+            pybamm.DiscretisationError, "Cannot automatically discretise model"
+        ):
+            solver.set_up(model, {})
+
     def test_convert_to_casadi_format(self):
         # Make sure model is converted to casadi format
         model = pybamm.BaseModel()