diff --git a/.gitignore b/.gitignore
index d59641db3..444611541 100644
--- a/.gitignore
+++ b/.gitignore
@@ -4,6 +4,7 @@
 *.png
 /local/
 *.DS_Store
+*.mat
 
 # don't ignore important .txt files
 !requirements*
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 9bcc47938..348939104 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -3,6 +3,7 @@
 ## Features
 
 -   Added `InputParameter` node for quickly changing parameter values ([#752](https://github.com/pybamm-team/PyBaMM/pull/752))
+-   Added optional R(x) distribution in particle models ([#745](https://github.com/pybamm-team/PyBaMM/pull/745))
 -   Generalized importing of external variables ([#728](https://github.com/pybamm-team/PyBaMM/pull/728))
 -   Separated active and inactive material volume fractions ([#726](https://github.com/pybamm-team/PyBaMM/pull/726))
 -   Added submodels for tortuosity ([#726](https://github.com/pybamm-team/PyBaMM/pull/726))
diff --git a/examples/scripts/SPMe_SOC.py b/examples/scripts/SPMe_SOC.py
index 411eb50f6..0ba2b5295 100644
--- a/examples/scripts/SPMe_SOC.py
+++ b/examples/scripts/SPMe_SOC.py
@@ -27,14 +27,14 @@ for l_n in thicknesses:
     w = 0.207 / factor
     A = h * w
     l_s = 2.5e-5
-    l1d = (l_n + l_p + l_s)
+    l1d = l_n + l_p + l_s
     vol = h * w * l1d
     vol_cm3 = vol * 1e6
     tot_cap = 0.0
     tot_time = 0.0
     fig, axes = plt.subplots(1, 2, sharey=True)
     I_mag = 1.01 / factor
-    print('*' * 30)
+    print("*" * 30)
     for enum, I_app in enumerate([-1.0, 1.0]):
         I_app *= I_mag
         # load model
@@ -44,27 +44,33 @@ for l_n in thicknesses:
         # load parameter values and process model and geometry
         param = model.default_parameter_values
         param.update(
-            {"Electrode height [m]": h,
-             "Electrode width [m]": w,
-             "Negative electrode thickness [m]": l_n,
-             "Positive electrode thickness [m]": l_p,
-             "Separator thickness [m]": l_s,
-             "Lower voltage cut-off [V]": 2.8,
-             "Upper voltage cut-off [V]": 4.7,
-             "Maximum concentration in negative electrode [mol.m-3]": 25000,
-             "Maximum concentration in positive electrode [mol.m-3]": 50000,
-             "Initial concentration in negative electrode [mol.m-3]": 12500,
-             "Initial concentration in positive electrode [mol.m-3]": 25000,
-             "Negative electrode surface area density [m-1]": 180000.0,
-             "Positive electrode surface area density [m-1]": 150000.0,
-             "Typical current [A]": I_app,
-             }
+            {
+                "Electrode height [m]": h,
+                "Electrode width [m]": w,
+                "Negative electrode thickness [m]": l_n,
+                "Positive electrode thickness [m]": l_p,
+                "Separator thickness [m]": l_s,
+                "Lower voltage cut-off [V]": 2.8,
+                "Upper voltage cut-off [V]": 4.7,
+                "Maximum concentration in negative electrode [mol.m-3]": 25000,
+                "Maximum concentration in positive electrode [mol.m-3]": 50000,
+                "Initial concentration in negative electrode [mol.m-3]": 12500,
+                "Initial concentration in positive electrode [mol.m-3]": 25000,
+                "Negative electrode surface area density [m-1]": 180000.0,
+                "Positive electrode surface area density [m-1]": 150000.0,
+                "Typical current [A]": I_app,
+            }
         )
         param.process_model(model)
         param.process_geometry(geometry)
         s_var = pybamm.standard_spatial_vars
-        var_pts = {s_var.x_n: 5, s_var.x_s: 5, s_var.x_p: 5,
-                   s_var.r_n: 5, s_var.r_p: 10}
+        var_pts = {
+            s_var.x_n: 5,
+            s_var.x_s: 5,
+            s_var.x_p: 5,
+            s_var.r_n: 5,
+            s_var.r_p: 10,
+        }
         # set mesh
         mesh = pybamm.Mesh(geometry, model.default_submesh_types, var_pts)
         # discretise model
@@ -74,58 +80,61 @@ for l_n in thicknesses:
         t_eval = np.linspace(0, 0.2, 100)
         sol = model.default_solver.solve(model, t_eval)
         var = "Positive electrode average extent of lithiation"
-        xpext = pybamm.ProcessedVariable(model.variables[var],
-                                         sol.t, sol.y, mesh=mesh)
+        xpext = pybamm.ProcessedVariable(model.variables[var], sol.t, sol.y, mesh=mesh)
         var = "Negative electrode average extent of lithiation"
-        xnext = pybamm.ProcessedVariable(model.variables[var],
-                                         sol.t, sol.y, mesh=mesh)
+        xnext = pybamm.ProcessedVariable(model.variables[var], sol.t, sol.y, mesh=mesh)
         var = "X-averaged positive particle surface concentration"
-        xpsurf = pybamm.ProcessedVariable(model.variables[var],
-                                          sol.t, sol.y, mesh=mesh)
+        xpsurf = pybamm.ProcessedVariable(model.variables[var], sol.t, sol.y, mesh=mesh)
         var = "X-averaged negative particle surface concentration"
-        xnsurf = pybamm.ProcessedVariable(model.variables[var],
-                                          sol.t, sol.y, mesh=mesh)
-        time = pybamm.ProcessedVariable(model.variables["Time [h]"],
-                                        sol.t, sol.y, mesh=mesh)
+        xnsurf = pybamm.ProcessedVariable(model.variables[var], sol.t, sol.y, mesh=mesh)
+        time = pybamm.ProcessedVariable(
+            model.variables["Time [h]"], sol.t, sol.y, mesh=mesh
+        )
         # Coulomb counting
         time_hours = time(sol.t)
         dc_time = np.around(time_hours[-1], 3)
         # Capacity mAh
         cap = np.absolute(I_app * 1000 * dc_time)
         cap_time = np.absolute(I_app * 1000 * time_hours)
-        axes[enum].plot(cap_time,
-                        xnext(sol.t), 'r-', label='Average Neg')
-        axes[enum].plot(cap_time,
-                        xpext(sol.t), 'b-', label='Average Pos')
-        axes[enum].plot(cap_time,
-                        xnsurf(sol.t), 'r--', label='Surface Neg')
-        axes[enum].plot(cap_time,
-                        xpsurf(sol.t), 'b--', label='Surface Pos')
-        axes[enum].set_xlabel('Capacity [mAh]')
+        axes[enum].plot(cap_time, xnext(sol.t), "r-", label="Average Neg")
+        axes[enum].plot(cap_time, xpext(sol.t), "b-", label="Average Pos")
+        axes[enum].plot(cap_time, xnsurf(sol.t), "r--", label="Surface Neg")
+        axes[enum].plot(cap_time, xpsurf(sol.t), "b--", label="Surface Pos")
+        axes[enum].set_xlabel("Capacity [mAh]")
         handles, labels = axes[enum].get_legend_handles_labels()
         axes[enum].legend(handles, labels)
         if I_app < 0.0:
-            axes[enum].set_ylabel('Extent of Lithiation, Elecrode Ratio: '
-                                  + str(e_ratio))
-            axes[enum].title.set_text('Charge')
+            axes[enum].set_ylabel(
+                "Extent of Lithiation, Elecrode Ratio: " + str(e_ratio)
+            )
+            axes[enum].title.set_text("Charge")
         else:
-            axes[enum].title.set_text('Discharge')
-        print('Applied Current', I_app, 'A', 'Time',
-              dc_time, 'hrs', 'Capacity', cap, 'mAh')
+            axes[enum].title.set_text("Discharge")
+        print(
+            "Applied Current",
+            I_app,
+            "A",
+            "Time",
+            dc_time,
+            "hrs",
+            "Capacity",
+            cap,
+            "mAh",
+        )
         tot_cap += cap
         tot_time += dc_time
 
-    print('Anode : Cathode thickness', e_ratio)
-    print('Total Charge/Discharge Time', tot_time, 'hrs')
-    print('Total Capacity', np.around(tot_cap, 3), 'mAh')
+    print("Anode : Cathode thickness", e_ratio)
+    print("Total Charge/Discharge Time", tot_time, "hrs")
+    print("Total Capacity", np.around(tot_cap, 3), "mAh")
     specific_cap = np.around(tot_cap, 3) / vol_cm3
-    print('Total Capacity', specific_cap, 'mAh.cm-3')
+    print("Total Capacity", specific_cap, "mAh.cm-3")
     capacities.append(tot_cap)
     specific_capacities.append(specific_cap)
 
 fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True)
 ax1.plot(thicknesses / l_p, capacities)
 ax2.plot(thicknesses / l_p, specific_capacities)
-ax1.set_ylabel('Capacity [mAh]')
-ax2.set_ylabel('Specific Capacity [mAh.cm-3]')
-ax2.set_xlabel('Anode : Cathode thickness')
+ax1.set_ylabel("Capacity [mAh]")
+ax2.set_ylabel("Specific Capacity [mAh.cm-3]")
+ax2.set_xlabel("Anode : Cathode thickness")
diff --git a/examples/scripts/SPMe_step.py b/examples/scripts/SPMe_step.py
index 666976c61..f669deefd 100644
--- a/examples/scripts/SPMe_step.py
+++ b/examples/scripts/SPMe_step.py
@@ -49,7 +49,6 @@ while time < end_time:
 # plot
 voltage = pybamm.ProcessedVariable(
     model.variables["Terminal voltage [V]"], solution.t, solution.y, mesh=mesh
-
 )
 step_voltage = pybamm.ProcessedVariable(
     model.variables["Terminal voltage [V]"], step_solution.t, step_solution.y, mesh=mesh
diff --git a/examples/scripts/compare_extrapolations.py b/examples/scripts/compare_extrapolations.py
index defcf557f..5ef9d9883 100644
--- a/examples/scripts/compare_extrapolations.py
+++ b/examples/scripts/compare_extrapolations.py
@@ -29,4 +29,3 @@ models = [sim_lin.built_model, sim_quad.built_model]
 solutions = [sim_lin.solution, sim_quad.solution]
 plot = pybamm.QuickPlot(models, sim_lin.mesh, solutions)
 plot.dynamic_plot()
-
diff --git a/examples/scripts/compare_lithium_ion_particle_distribution.py b/examples/scripts/compare_lithium_ion_particle_distribution.py
new file mode 100644
index 000000000..57d68d089
--- /dev/null
+++ b/examples/scripts/compare_lithium_ion_particle_distribution.py
@@ -0,0 +1,82 @@
+#
+# Compare lithium-ion battery models
+#
+import argparse
+import numpy as np
+import pybamm
+
+parser = argparse.ArgumentParser()
+parser.add_argument(
+    "--debug", action="store_true", help="Set logging level to 'DEBUG'."
+)
+args = parser.parse_args()
+if args.debug:
+    pybamm.set_logging_level("DEBUG")
+else:
+    pybamm.set_logging_level("INFO")
+
+# load models
+options = {"thermal": "isothermal"}
+models = [
+    pybamm.lithium_ion.DFN(options, name="standard DFN"),
+    pybamm.lithium_ion.DFN(options, name="particle DFN"),
+]
+
+
+# load parameter values and process models and geometry
+params = [models[0].default_parameter_values, models[1].default_parameter_values]
+params[0]["Typical current [A]"] = 1.0
+params[0].process_model(models[0])
+
+
+params[1]["Typical current [A]"] = 1.0
+
+
+def negative_distribution(x):
+    return 1 + x
+
+
+def positive_distribution(x):
+    return 1 + (x - (1 - models[1].param.l_p))
+
+
+params[1]["Negative particle distribution in x"] = negative_distribution
+params[1]["Positive particle distribution in x"] = positive_distribution
+params[1].process_model(models[1])
+
+# set mesh
+var = pybamm.standard_spatial_vars
+var_pts = {var.x_n: 10, var.x_s: 10, var.x_p: 10, var.r_n: 5, var.r_p: 5}
+
+# discretise models
+for param, model in zip(params, models):
+    # create geometry
+    geometry = model.default_geometry
+    param.process_geometry(geometry)
+    mesh = pybamm.Mesh(geometry, models[-1].default_submesh_types, var_pts)
+    disc = pybamm.Discretisation(mesh, model.default_spatial_methods)
+    disc.process_model(model)
+
+# solve model
+solutions = [None] * len(models)
+t_eval = np.linspace(0, 0.3, 100)
+for i, model in enumerate(models):
+    solutions[i] = model.default_solver.solve(model, t_eval)
+
+
+output_variables = [
+    "Negative particle surface concentration",
+    "Electrolyte concentration",
+    "Positive particle surface concentration",
+    "Current [A]",
+    "Negative electrode potential [V]",
+    "Electrolyte potential [V]",
+    "Positive electrode potential [V]",
+    "Terminal voltage [V]",
+    "Negative particle distribution in x",
+    "Positive particle distribution in x",
+]
+
+# plot
+plot = pybamm.QuickPlot(models, mesh, solutions, output_variables=output_variables)
+plot.dynamic_plot()
diff --git a/examples/scripts/heat_equation.py b/examples/scripts/heat_equation.py
index 5bccd4fd9..8ea4a2a0c 100644
--- a/examples/scripts/heat_equation.py
+++ b/examples/scripts/heat_equation.py
@@ -122,11 +122,7 @@ for i, t in enumerate(plot_times):
         label="Numerical" if i == 0 else "",
     )
     plt.plot(
-        xx,
-        T_exact(xx, t),
-        "-",
-        color=color,
-        label="Exact (t={})".format(plot_times[i]),
+        xx, T_exact(xx, t), "-", color=color, label="Exact (t={})".format(plot_times[i])
     )
 plt.xlabel("x", fontsize=16)
 plt.ylabel("T", fontsize=16)
diff --git a/input/parameters/lithium-ion/anodes/graphite_mcmb2528_Marquis2019/parameters.csv b/input/parameters/lithium-ion/anodes/graphite_mcmb2528_Marquis2019/parameters.csv
index 3d1469d95..9a5fcbcd2 100644
--- a/input/parameters/lithium-ion/anodes/graphite_mcmb2528_Marquis2019/parameters.csv
+++ b/input/parameters/lithium-ion/anodes/graphite_mcmb2528_Marquis2019/parameters.csv
@@ -11,6 +11,7 @@ Negative electrode OCP [V],[function]graphite_mcmb2528_ocp_Dualfoil1998,
 Negative electrode porosity,0.3,Scott Moura FastDFN,electrolyte volume fraction
 Negative electrode active material volume fraction,0.7,,assuming zero binder volume fraction
 Negative particle radius [m],1E-05,Scott Moura FastDFN,
+Negative particle distribution in x,1,,
 Negative electrode surface area density [m-1],180000,Scott Moura FastDFN,
 Negative electrode Bruggeman coefficient (electrolyte),1.5,Scott Moura FastDFN,
 Negative electrode Bruggeman coefficient (electrode),1.5,Scott Moura FastDFN,
diff --git a/input/parameters/lithium-ion/cathodes/lico2_Marquis2019/parameters.csv b/input/parameters/lithium-ion/cathodes/lico2_Marquis2019/parameters.csv
index d1fabddca..1eafb9b65 100644
--- a/input/parameters/lithium-ion/cathodes/lico2_Marquis2019/parameters.csv
+++ b/input/parameters/lithium-ion/cathodes/lico2_Marquis2019/parameters.csv
@@ -11,6 +11,7 @@ Positive electrode OCP [V],[function]lico2_ocp_Dualfoil1998,
 Positive electrode porosity,0.3,Scott Moura FastDFN,electrolyte volume fraction
 Positive electrode active material volume fraction,0.7,,assuming zero binder volume fraction
 Positive particle radius [m],1E-05,Scott Moura FastDFN,
+Positive particle distribution in x,1,,
 Positive electrode surface area density [m-1],150000,Scott Moura FastDFN,
 Positive electrode Bruggeman coefficient (electrolyte),1.5,Scott Moura FastDFN,
 Positive electrode Bruggeman coefficient (electrode),1.5,Scott Moura FastDFN,
diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index 5c5a18353..a97a49f57 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -100,7 +100,12 @@ from .expression_tree.functions import *
 from .expression_tree.interpolant import Interpolant
 from .expression_tree.input_parameter import InputParameter
 from .expression_tree.parameter import Parameter, FunctionParameter
-from .expression_tree.broadcasts import Broadcast, PrimaryBroadcast, FullBroadcast
+from .expression_tree.broadcasts import (
+    Broadcast,
+    PrimaryBroadcast,
+    FullBroadcast,
+    ones_like,
+)
 from .expression_tree.scalar import Scalar
 from .expression_tree.variable import Variable
 from .expression_tree.independent_variable import (
diff --git a/pybamm/expression_tree/broadcasts.py b/pybamm/expression_tree/broadcasts.py
index 5395b880c..5e32b2809 100644
--- a/pybamm/expression_tree/broadcasts.py
+++ b/pybamm/expression_tree/broadcasts.py
@@ -184,3 +184,26 @@ class FullBroadcast(Broadcast):
         )
 
         return child_eval * vec
+
+
+def ones_like(*symbols):
+    """
+    Create a symbol with the same shape as the input symbol and with constant value '1',
+    using `FullBroadcast`.
+
+    Parameters
+    ----------
+    symbols : :class:`Symbol`
+        Symbols whose shape to copy
+    """
+    # Make a symbol that combines all the children, to get the right domain
+    # that takes all the child symbols into account
+    sum_symbol = symbols[0]
+    for sym in symbols:
+        sum_symbol += sym
+
+    # Just return scalar 1 if symbol has no domain (no broadcasting necessary)
+    if sum_symbol.domain == []:
+        return pybamm.Scalar(1)
+    else:
+        return FullBroadcast(1, sum_symbol.domain, sum_symbol.auxiliary_domains)
diff --git a/pybamm/models/submodels/particle/fickian/base_fickian_particle.py b/pybamm/models/submodels/particle/fickian/base_fickian_particle.py
index 216164a63..368927c5c 100644
--- a/pybamm/models/submodels/particle/fickian/base_fickian_particle.py
+++ b/pybamm/models/submodels/particle/fickian/base_fickian_particle.py
@@ -35,16 +35,6 @@ class BaseModel(BaseParticle):
     def _unpack(self, variables):
         raise NotImplementedError
 
-    def set_rhs(self, variables):
-
-        c, N, _ = self._unpack(variables)
-
-        if self.domain == "Negative":
-            self.rhs = {c: -(1 / self.param.C_n) * pybamm.div(N)}
-
-        elif self.domain == "Positive":
-            self.rhs = {c: -(1 / self.param.C_p) * pybamm.div(N)}
-
     def set_boundary_conditions(self, variables):
 
         c, _, j = self._unpack(variables)
diff --git a/pybamm/models/submodels/particle/fickian/fickian_many_particles.py b/pybamm/models/submodels/particle/fickian/fickian_many_particles.py
index 159a0a710..8b0d9d40c 100644
--- a/pybamm/models/submodels/particle/fickian/fickian_many_particles.py
+++ b/pybamm/models/submodels/particle/fickian/fickian_many_particles.py
@@ -46,8 +46,30 @@ class ManyParticles(BaseModel):
         N_s = self._flux_law(c_s, T_k)
 
         variables.update(self._get_standard_flux_variables(N_s, N_s))
+
+        if self.domain == "Negative":
+            x = pybamm.standard_spatial_vars.x_n
+            R = pybamm.FunctionParameter("Negative particle distribution in x", x)
+            variables.update({"Negative particle distribution in x": R})
+
+        elif self.domain == "Positive":
+            x = pybamm.standard_spatial_vars.x_p
+            R = pybamm.FunctionParameter("Positive particle distribution in x", x)
+            variables.update({"Positive particle distribution in x": R})
         return variables
 
+    def set_rhs(self, variables):
+
+        c, N, _ = self._unpack(variables)
+
+        if self.domain == "Negative":
+            R = variables["Negative particle distribution in x"]
+            self.rhs = {c: -(1 / (R ** 2 * self.param.C_n)) * pybamm.div(N)}
+
+        elif self.domain == "Positive":
+            R = variables["Positive particle distribution in x"]
+            self.rhs = {c: -(1 / (R ** 2 * self.param.C_p)) * pybamm.div(N)}
+
     def _unpack(self, variables):
         c_s = variables[self.domain + " particle concentration"]
         N_s = variables[self.domain + " particle flux"]
diff --git a/pybamm/models/submodels/particle/fickian/fickian_single_particle.py b/pybamm/models/submodels/particle/fickian/fickian_single_particle.py
index db3486a9d..51bed8f20 100644
--- a/pybamm/models/submodels/particle/fickian/fickian_single_particle.py
+++ b/pybamm/models/submodels/particle/fickian/fickian_single_particle.py
@@ -52,6 +52,16 @@ class SingleParticle(BaseModel):
 
         return variables
 
+    def set_rhs(self, variables):
+
+        c, N, _ = self._unpack(variables)
+
+        if self.domain == "Negative":
+            self.rhs = {c: -(1 / self.param.C_n) * pybamm.div(N)}
+
+        elif self.domain == "Positive":
+            self.rhs = {c: -(1 / self.param.C_p) * pybamm.div(N)}
+
     def _unpack(self, variables):
         c_s_xav = variables[
             "X-averaged " + self.domain.lower() + " particle concentration"
diff --git a/pybamm/parameters/parameter_values.py b/pybamm/parameters/parameter_values.py
index 2e98c7270..488b68f28 100644
--- a/pybamm/parameters/parameter_values.py
+++ b/pybamm/parameters/parameter_values.py
@@ -441,8 +441,11 @@ class ParameterValues(dict):
                 function = pybamm.Interpolant(data, *new_children, name=name)
             elif isinstance(function_name, numbers.Number):
                 # If the "function" is provided is actually a scalar, return a Scalar
-                # object instead of throwing an error
-                function = pybamm.Scalar(function_name, name=symbol.name)
+                # object instead of throwing an error.
+                # Also use ones_like so that we get the right shapes
+                function = pybamm.Scalar(
+                    function_name, name=symbol.name
+                ) * pybamm.ones_like(*new_children)
             else:
                 # otherwise evaluate the function to create a new PyBaMM object
                 function = function_name(*new_children)

Particle Size Distribution for DFN
**Summary**
Adding particle size distribution in x for DFN model

 **Motivation**
Particle size is actually different from particle to particle.

**Additional context**

diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
index b8ac65125..938df0a3d 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py
@@ -95,6 +95,21 @@ class TestDFN(unittest.TestCase):
         modeltest = tests.StandardModelTest(model)
         modeltest.test_all()
 
+    def test_particle_distribution_in_x(self):
+        model = pybamm.lithium_ion.DFN()
+        param = model.default_parameter_values
+
+        def negative_distribution(x):
+            return 1 + x
+
+        def positive_distribution(x):
+            return 1 + (x - (1 - model.param.l_p))
+
+        param["Negative particle distribution in x"] = negative_distribution
+        param["Positive particle distribution in x"] = positive_distribution
+        modeltest = tests.StandardModelTest(model, parameter_values=param)
+        modeltest.test_all()
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/unit/test_expression_tree/test_broadcasts.py b/tests/unit/test_expression_tree/test_broadcasts.py
index a331652d9..c91a45386 100644
--- a/tests/unit/test_expression_tree/test_broadcasts.py
+++ b/tests/unit/test_expression_tree/test_broadcasts.py
@@ -26,6 +26,29 @@ class TestBroadcasts(unittest.TestCase):
         with self.assertRaisesRegex(ValueError, "Variables on the current collector"):
             pybamm.Broadcast(a, "electrode")
 
+    def test_ones_like(self):
+        a = pybamm.Variable("a")
+        ones_like_a = pybamm.ones_like(a)
+        self.assertEqual(ones_like_a.id, pybamm.Scalar(1).id)
+
+        a = pybamm.Variable(
+            "a",
+            domain="negative electrode",
+            auxiliary_domains={"secondary": "current collector"},
+        )
+        ones_like_a = pybamm.ones_like(a)
+        self.assertIsInstance(ones_like_a, pybamm.FullBroadcast)
+        self.assertEqual(ones_like_a.name, "broadcast")
+        self.assertEqual(ones_like_a.domain, a.domain)
+        self.assertEqual(ones_like_a.auxiliary_domains, a.auxiliary_domains)
+
+        b = pybamm.Variable("b", domain="current collector")
+        ones_like_ab = pybamm.ones_like(b, a)
+        self.assertIsInstance(ones_like_ab, pybamm.FullBroadcast)
+        self.assertEqual(ones_like_ab.name, "broadcast")
+        self.assertEqual(ones_like_ab.domain, a.domain)
+        self.assertEqual(ones_like_ab.auxiliary_domains, a.auxiliary_domains)
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/unit/test_parameters/test_parameter_values.py b/tests/unit/test_parameters/test_parameter_values.py
index c5182c3c8..225303211 100644
--- a/tests/unit/test_parameters/test_parameter_values.py
+++ b/tests/unit/test_parameters/test_parameter_values.py
@@ -263,7 +263,9 @@ class TestParameterValues(unittest.TestCase):
         # process constant function
         const = pybamm.FunctionParameter("const", a)
         processed_const = parameter_values.process_symbol(const)
-        self.assertIsInstance(processed_const, pybamm.Scalar)
+        self.assertIsInstance(processed_const, pybamm.Multiplication)
+        self.assertIsInstance(processed_const.left, pybamm.Scalar)
+        self.assertIsInstance(processed_const.right, pybamm.Scalar)
         self.assertEqual(processed_const.evaluate(), 254)
 
         # process differentiated function parameter
diff --git a/tests/unit/test_parameters/test_update_parameters.py b/tests/unit/test_parameters/test_update_parameters.py
index e26c153da..b3e80b790 100644
--- a/tests/unit/test_parameters/test_update_parameters.py
+++ b/tests/unit/test_parameters/test_update_parameters.py
@@ -73,7 +73,9 @@ class TestUpdateParameters(unittest.TestCase):
         modeltest3.test_solving(t_eval=t_eval)
         Y3 = modeltest3.solution.y
 
-        self.assertIsInstance(model3.variables["Current [A]"], pybamm.Scalar)
+        self.assertIsInstance(model3.variables["Current [A]"], pybamm.Multiplication)
+        self.assertIsInstance(model3.variables["Current [A]"].left, pybamm.Scalar)
+        self.assertIsInstance(model3.variables["Current [A]"].right, pybamm.Scalar)
         self.assertEqual(model3.variables["Current [A]"].evaluate(), 0.0)
 
         # results should be different

pip install -e .[dev]
parse_log_pytest
requirements.txt
pytest
apt-get update apt-get install -y gcc
3.9
requirements/base.txt
pytest --no-header -rA --tb=line --color=no -p no:cacheprovider -W ignore::DeprecationWarning
anyio==4.9.0
anytree==2.12.1
argon2-cffi==23.1.0
argon2-cffi-bindings==21.2.0
arrow==1.3.0
asttokens==3.0.0
async-lru==2.0.5
attrs==25.3.0
autograd==1.7.0
babel==2.17.0
beautifulsoup4==4.13.3
black==25.1.0
bleach==6.2.0
casadi==3.7.0
certifi==2025.1.31
cffi==1.17.1
charset-normalizer==3.4.1
click==8.1.8
comm==0.2.2
contourpy==1.3.0
cycler==0.12.1
debugpy==1.8.13
decorator==5.2.1
defusedxml==0.7.1
exceptiongroup==1.2.2
executing==2.2.0
fastjsonschema==2.21.1
flake8==7.2.0
fonttools==4.56.0
fqdn==1.5.1
h11==0.14.0
httpcore==1.0.7
httpx==0.28.1
idna==3.10
importlib_metadata==8.6.1
importlib_resources==6.5.2
iniconfig==2.1.0
ipykernel==6.29.5
ipython==8.18.1
ipywidgets==8.1.5
isoduration==20.11.0
jedi==0.19.2
Jinja2==3.1.6
json5==0.10.0
jsonpointer==3.0.0
jsonschema==4.23.0
jsonschema-specifications==2024.10.1
jupyter==1.1.1
jupyter-console==6.6.3
jupyter-events==0.12.0
jupyter-lsp==2.2.5
jupyter_client==8.6.3
jupyter_core==5.7.2
jupyter_server==2.15.0
jupyter_server_terminals==0.5.3
jupyterlab==4.3.6
jupyterlab_pygments==0.3.0
jupyterlab_server==2.27.3
jupyterlab_widgets==3.0.13
kiwisolver==1.4.7
MarkupSafe==3.0.2
matplotlib==3.9.4
matplotlib-inline==0.1.7
mccabe==0.7.0
mistune==3.1.3
mypy-extensions==1.0.0
nbclient==0.10.2
nbconvert==7.16.6
nbformat==5.10.4
nest-asyncio==1.6.0
notebook==7.3.3
notebook_shim==0.2.4
numpy==2.0.2
overrides==7.7.0
packaging==24.2
pandas==2.2.3
pandocfilters==1.5.1
parso==0.8.4
pathspec==0.12.1
pexpect==4.9.0
pillow==11.1.0
platformdirs==4.3.7
pluggy==1.5.0
prometheus_client==0.21.1
prompt_toolkit==3.0.50
psutil==7.0.0
ptyprocess==0.7.0
pure_eval==0.2.3
-e git+https://github.com/pybamm-team/PyBaMM.git@61561b0ed9794364723ce1a1f150eb222214603d#egg=pybamm
pycodestyle==2.13.0
pycparser==2.22
pyflakes==3.3.1
Pygments==2.19.1
pyparsing==3.2.3
pytest==8.3.5
python-dateutil==2.9.0.post0
python-json-logger==3.3.0
pytz==2025.2
PyYAML==6.0.2
pyzmq==26.3.0
referencing==0.36.2
requests==2.32.3
rfc3339-validator==0.1.4
rfc3986-validator==0.1.1
rpds-py==0.24.0
scikit-fem==10.0.2
scipy==1.13.1
Send2Trash==1.8.3
six==1.17.0
sniffio==1.3.1
soupsieve==2.6
stack-data==0.6.3
terminado==0.18.1
tinycss2==1.4.0
tomli==2.2.1
tornado==6.4.2
traitlets==5.14.3
types-python-dateutil==2.9.0.20241206
typing_extensions==4.13.0
tzdata==2025.2
uri-template==1.3.0
urllib3==2.3.0
wcwidth==0.2.13
webcolors==24.11.1
webencodings==0.5.1
websocket-client==1.8.0
widgetsnbextension==4.0.13
zipp==3.21.0

name: PyBaMM
channels:
  - defaults
  - https://repo.anaconda.com/pkgs/main
  - https://repo.anaconda.com/pkgs/r
  - conda-forge
dependencies:
  - _libgcc_mutex=0.1=main
  - _openmp_mutex=5.1=1_gnu
  - ca-certificates=2025.2.25=h06a4308_0
  - ld_impl_linux-64=2.40=h12ee557_0
  - libffi=3.4.4=h6a678d5_1
  - libgcc-ng=11.2.0=h1234567_1
  - libgomp=11.2.0=h1234567_1
  - libstdcxx-ng=11.2.0=h1234567_1
  - ncurses=6.4=h6a678d5_0
  - openssl=3.0.16=h5eee18b_0
  - pip=25.0=py39h06a4308_0
  - python=3.9.21=he870216_1
  - readline=8.2=h5eee18b_0
  - setuptools=75.8.0=py39h06a4308_0
  - sqlite=3.45.3=h5eee18b_0
  - tk=8.6.14=h39e8969_0
  - wheel=0.45.1=py39h06a4308_0
  - xz=5.6.4=h5eee18b_1
  - zlib=1.2.13=h5eee18b_1
  - pip:
      - anyio==4.9.0
      - anytree==2.12.1
      - argon2-cffi==23.1.0
      - argon2-cffi-bindings==21.2.0
      - arrow==1.3.0
      - asttokens==3.0.0
      - async-lru==2.0.5
      - attrs==25.3.0
      - autograd==1.7.0
      - babel==2.17.0
      - beautifulsoup4==4.13.3
      - black==25.1.0
      - bleach==6.2.0
      - casadi==3.7.0
      - certifi==2025.1.31
      - cffi==1.17.1
      - charset-normalizer==3.4.1
      - click==8.1.8
      - comm==0.2.2
      - contourpy==1.3.0
      - cycler==0.12.1
      - debugpy==1.8.13
      - decorator==5.2.1
      - defusedxml==0.7.1
      - exceptiongroup==1.2.2
      - executing==2.2.0
      - fastjsonschema==2.21.1
      - flake8==7.2.0
      - fonttools==4.56.0
      - fqdn==1.5.1
      - h11==0.14.0
      - httpcore==1.0.7
      - httpx==0.28.1
      - idna==3.10
      - importlib-metadata==8.6.1
      - importlib-resources==6.5.2
      - iniconfig==2.1.0
      - ipykernel==6.29.5
      - ipython==8.18.1
      - ipywidgets==8.1.5
      - isoduration==20.11.0
      - jedi==0.19.2
      - jinja2==3.1.6
      - json5==0.10.0
      - jsonpointer==3.0.0
      - jsonschema==4.23.0
      - jsonschema-specifications==2024.10.1
      - jupyter==1.1.1
      - jupyter-client==8.6.3
      - jupyter-console==6.6.3
      - jupyter-core==5.7.2
      - jupyter-events==0.12.0
      - jupyter-lsp==2.2.5
      - jupyter-server==2.15.0
      - jupyter-server-terminals==0.5.3
      - jupyterlab==4.3.6
      - jupyterlab-pygments==0.3.0
      - jupyterlab-server==2.27.3
      - jupyterlab-widgets==3.0.13
      - kiwisolver==1.4.7
      - markupsafe==3.0.2
      - matplotlib==3.9.4
      - matplotlib-inline==0.1.7
      - mccabe==0.7.0
      - mistune==3.1.3
      - mypy-extensions==1.0.0
      - nbclient==0.10.2
      - nbconvert==7.16.6
      - nbformat==5.10.4
      - nest-asyncio==1.6.0
      - notebook==7.3.3
      - notebook-shim==0.2.4
      - numpy==2.0.2
      - overrides==7.7.0
      - packaging==24.2
      - pandas==2.2.3
      - pandocfilters==1.5.1
      - parso==0.8.4
      - pathspec==0.12.1
      - pexpect==4.9.0
      - pillow==11.1.0
      - platformdirs==4.3.7
      - pluggy==1.5.0
      - prometheus-client==0.21.1
      - prompt-toolkit==3.0.50
      - psutil==7.0.0
      - ptyprocess==0.7.0
      - pure-eval==0.2.3
      - pycodestyle==2.13.0
      - pycparser==2.22
      - pyflakes==3.3.1
      - pygments==2.19.1
      - pyparsing==3.2.3
      - pytest==8.3.5
      - python-dateutil==2.9.0.post0
      - python-json-logger==3.3.0
      - pytz==2025.2
      - pyyaml==6.0.2
      - pyzmq==26.3.0
      - referencing==0.36.2
      - requests==2.32.3
      - rfc3339-validator==0.1.4
      - rfc3986-validator==0.1.1
      - rpds-py==0.24.0
      - scikit-fem==10.0.2
      - scipy==1.13.1
      - send2trash==1.8.3
      - six==1.17.0
      - sniffio==1.3.1
      - soupsieve==2.6
      - stack-data==0.6.3
      - terminado==0.18.1
      - tinycss2==1.4.0
      - tomli==2.2.1
      - tornado==6.4.2
      - traitlets==5.14.3
      - types-python-dateutil==2.9.0.20241206
      - typing-extensions==4.13.0
      - tzdata==2025.2
      - uri-template==1.3.0
      - urllib3==2.3.0
      - wcwidth==0.2.13
      - webcolors==24.11.1
      - webencodings==0.5.1
      - websocket-client==1.8.0
      - widgetsnbextension==4.0.13
      - zipp==3.21.0
prefix: /opt/conda/envs/PyBaMM

tests/unit/test_expression_tree/test_broadcasts.py::TestBroadcasts::test_ones_like tests/unit/test_parameters/test_parameter_values.py::TestParameterValues::test_process_function_parameter tests/unit/test_parameters/test_update_parameters.py::TestUpdateParameters::test_update_model
tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py::TestDFN::test_basic_processing_2plus1D
tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py::TestDFN::test_basic_processing tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py::TestDFN::test_basic_processing_1plus1D tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py::TestDFN::test_full_thermal tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py::TestDFN::test_lumped_thermal tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py::TestDFN::test_optimisations tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py::TestDFN::test_particle_distribution_in_x tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py::TestDFN::test_particle_fast_diffusion tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py::TestDFN::test_set_up tests/unit/test_expression_tree/test_broadcasts.py::TestBroadcasts::test_broadcast tests/unit/test_expression_tree/test_broadcasts.py::TestBroadcasts::test_broadcast_number tests/unit/test_expression_tree/test_broadcasts.py::TestBroadcasts::test_broadcast_type tests/unit/test_parameters/test_parameter_values.py::TestParameterValues::test_check_and_update_parameter_values tests/unit/test_parameters/test_parameter_values.py::TestParameterValues::test_evaluate tests/unit/test_parameters/test_parameter_values.py::TestParameterValues::test_init tests/unit/test_parameters/test_parameter_values.py::TestParameterValues::test_interpolant_against_function tests/unit/test_parameters/test_parameter_values.py::TestParameterValues::test_multi_var_function_parameter tests/unit/test_parameters/test_parameter_values.py::TestParameterValues::test_multi_var_function_with_parameters tests/unit/test_parameters/test_parameter_values.py::TestParameterValues::test_process_complex_expression tests/unit/test_parameters/test_parameter_values.py::TestParameterValues::test_process_empty_model tests/unit/test_parameters/test_parameter_values.py::TestParameterValues::test_process_inline_function_parameters tests/unit/test_parameters/test_parameter_values.py::TestParameterValues::test_process_input_parameter tests/unit/test_parameters/test_parameter_values.py::TestParameterValues::test_process_interpolant tests/unit/test_parameters/test_parameter_values.py::TestParameterValues::test_process_model tests/unit/test_parameters/test_parameter_values.py::TestParameterValues::test_process_symbol tests/unit/test_parameters/test_parameter_values.py::TestParameterValues::test_read_parameters_csv tests/unit/test_parameters/test_parameter_values.py::TestParameterValues::test_update tests/unit/test_parameters/test_parameter_values.py::TestParameterValues::test_update_from_chemistry tests/unit/test_parameters/test_update_parameters.py::TestUpdateParameters::test_inplace tests/unit/test_parameters/test_update_parameters.py::TestUpdateParameters::test_set_and_update_parameters tests/unit/test_parameters/test_update_parameters.py::TestUpdateParameters::test_update_geometry tests/unit/test_parameters/test_update_parameters.py::TestUpdateParameters::test_update_parameters_eqn
input/parameters/lithium-ion/anodes/graphite_mcmb2528_Marquis2019/parameters.csv pybamm/__init__.py pybamm/parameters/parameter_values.py pybamm/expression_tree/broadcasts.py examples/scripts/SPMe_step.py pybamm/models/submodels/particle/fickian/fickian_many_particles.py CHANGELOG.md pybamm/models/submodels/particle/fickian/fickian_single_particle.py .gitignore pybamm/models/submodels/particle/fickian/base_fickian_particle.py input/parameters/lithium-ion/cathodes/lico2_Marquis2019/parameters.csv examples/scripts/compare_extrapolations.py examples/scripts/compare_lithium_ion_particle_distribution.py examples/scripts/SPMe_SOC.py examples/scripts/heat_equation.py
input/parameters/lithium-ion/anodes/graphite_mcmb2528_Marquis2019/parameters.csv pybamm/__init__.py pybamm/parameters/parameter_values.py pybamm/expression_tree/broadcasts.py examples/scripts/SPMe_step.py pybamm/models/submodels/particle/fickian/fickian_many_particles.py CHANGELOG.md pybamm/models/submodels/particle/fickian/fickian_single_particle.py .gitignore pybamm/models/submodels/particle/fickian/base_fickian_particle.py input/parameters/lithium-ion/cathodes/lico2_Marquis2019/parameters.csv examples/scripts/compare_extrapolations.py examples/scripts/compare_lithium_ion_particle_distribution.py examples/scripts/SPMe_SOC.py examples/scripts/heat_equation.py