diff --git a/CHANGELOG.md b/CHANGELOG.md
index 4ffe632a8..a0c480ad9 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,6 +4,7 @@
 
 -   Reformatted Landesfeind electrolytes ([#1064](https://github.com/pybamm-team/PyBaMM/pull/1064))
 -   Adapted examples to be run in Google Colab ([#1061](https://github.com/pybamm-team/PyBaMM/pull/1061))
+-   Added some new solvers for algebraic models ([#1059](https://github.com/pybamm-team/PyBaMM/pull/1059))
 -   Added `length_scales` attribute to models ([#1058](https://github.com/pybamm-team/PyBaMM/pull/1058))
 -   Added averaging in secondary dimensions ([#1057](https://github.com/pybamm-team/PyBaMM/pull/1057))
 
diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index 2afdb9102..1e6de70f9 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -160,8 +160,6 @@ class Discretisation(object):
         # Set the y split for variables
         pybamm.logger.info("Set variable slices for {}".format(model.name))
         self.set_variable_slices(variables)
-        # Keep a record of y_slices in the model
-        model.y_slices = self.y_slices_explicit
 
         # now add extrapolated external variables to the boundary conditions
         # if required by the spatial method
@@ -183,6 +181,11 @@ class Discretisation(object):
             # create an empty copy of the original model
             model_disc = model.new_copy()
 
+        # Keep a record of y_slices in the model
+        model_disc.y_slices = self.y_slices_explicit
+        # Keep a record of the bounds in the model
+        model_disc.bounds = self.bounds
+
         model_disc.bcs = self.bcs
 
         pybamm.logger.info("Discretise initial conditions for {}".format(model.name))
@@ -240,11 +243,13 @@ class Discretisation(object):
         variables : iterable of :class:`pybamm.Variables`
             The variables for which to set slices
         """
-        # Set up y_slices
+        # Set up y_slices and bounds
         y_slices = defaultdict(list)
         y_slices_explicit = defaultdict(list)
         start = 0
         end = 0
+        lower_bounds = []
+        upper_bounds = []
         # Iterate through unpacked variables, adding appropriate slices to y_slices
         for variable in variables:
             # Add up the size of all the domains in variable.domain
@@ -261,13 +266,23 @@ class Discretisation(object):
                     for child, mesh in meshes.items():
                         for domain_mesh in mesh:
                             end += domain_mesh.npts_for_broadcast_to_nodes
+                        # Add to slices
                         y_slices[child.id].append(slice(start, end))
                         y_slices_explicit[child].append(slice(start, end))
+                        # Add to bounds
+                        lower_bounds.extend([child.bounds[0]] * (end - start))
+                        upper_bounds.extend([child.bounds[1]] * (end - start))
+                        # Increment start
                         start = end
             else:
                 end += self._get_variable_size(variable)
+                # Add to slices
                 y_slices[variable.id].append(slice(start, end))
                 y_slices_explicit[variable].append(slice(start, end))
+                # Add to bounds
+                lower_bounds.extend([variable.bounds[0]] * (end - start))
+                upper_bounds.extend([variable.bounds[1]] * (end - start))
+                # Increment start
                 start = end
 
         # Convert y_slices back to normal dictionary
@@ -275,6 +290,9 @@ class Discretisation(object):
         # Also keep a record of what the y_slices are, to be stored in the model
         self.y_slices_explicit = dict(y_slices_explicit)
 
+        # Also keep a record of bounds
+        self.bounds = (np.array(lower_bounds), np.array(upper_bounds))
+
         # reset discretised_symbols
         self._discretised_symbols = {}
 
@@ -885,7 +903,7 @@ class Discretisation(object):
                 # Broadcast new_child to the domain specified by symbol.domain
                 # Different discretisations may broadcast differently
                 if symbol.domain == []:
-                    symbol = disc_child * pybamm.Vector(np.array([1]))
+                    symbol = disc_child * pybamm.Vector([1])
                 else:
                     symbol = spatial_method.broadcast(
                         disc_child,
@@ -921,7 +939,7 @@ class Discretisation(object):
             return pybamm.StateVectorDot(
                 *self.y_slices[symbol.get_variable().id],
                 domain=symbol.domain,
-                auxiliary_domains=symbol.auxiliary_domains
+                auxiliary_domains=symbol.auxiliary_domains,
             )
 
         elif isinstance(symbol, pybamm.Variable):
@@ -972,7 +990,7 @@ class Discretisation(object):
                 return pybamm.StateVector(
                     *y_slices,
                     domain=symbol.domain,
-                    auxiliary_domains=symbol.auxiliary_domains
+                    auxiliary_domains=symbol.auxiliary_domains,
                 )
 
         elif isinstance(symbol, pybamm.SpatialVariable):
diff --git a/pybamm/expression_tree/array.py b/pybamm/expression_tree/array.py
index 41bed0708..95622571b 100644
--- a/pybamm/expression_tree/array.py
+++ b/pybamm/expression_tree/array.py
@@ -13,8 +13,9 @@ class Array(pybamm.Symbol):
     Parameters
     ----------
 
-    entries : numpy.array
-        the array associated with the node
+    entries : numpy.array or list
+        the array associated with the node. If a list is provided, it is converted to a
+        numpy array
     name : str, optional
         the name of the node
     domain : iterable of str, optional
@@ -35,6 +36,8 @@ class Array(pybamm.Symbol):
         auxiliary_domains=None,
         entries_string=None,
     ):
+        if isinstance(entries, list):
+            entries = np.array(entries)
         if entries.ndim == 1:
             entries = entries[:, np.newaxis]
         if name is None:
diff --git a/pybamm/expression_tree/concatenations.py b/pybamm/expression_tree/concatenations.py
index 611b45773..456029314 100644
--- a/pybamm/expression_tree/concatenations.py
+++ b/pybamm/expression_tree/concatenations.py
@@ -127,7 +127,7 @@ class NumpyConcatenation(Concatenation):
         # so that we can concatenate them
         for i, child in enumerate(children):
             if child.evaluates_to_number():
-                children[i] = child * pybamm.Vector(np.array([1]))
+                children[i] = child * pybamm.Vector([1])
         super().__init__(
             *children,
             name="numpy concatenation",
diff --git a/pybamm/expression_tree/matrix.py b/pybamm/expression_tree/matrix.py
index cd613266d..8aed40bc1 100644
--- a/pybamm/expression_tree/matrix.py
+++ b/pybamm/expression_tree/matrix.py
@@ -2,7 +2,7 @@
 # Matrix class
 #
 import pybamm
-
+import numpy as np
 from scipy.sparse import issparse
 
 
@@ -21,6 +21,8 @@ class Matrix(pybamm.Array):
         auxiliary_domains=None,
         entries_string=None,
     ):
+        if isinstance(entries, list):
+            entries = np.array(entries)
         if name is None:
             name = "Matrix {!s}".format(entries.shape)
             if issparse(entries):
diff --git a/pybamm/expression_tree/state_vector.py b/pybamm/expression_tree/state_vector.py
index 47a6c28e1..7bca90e6c 100644
--- a/pybamm/expression_tree/state_vector.py
+++ b/pybamm/expression_tree/state_vector.py
@@ -26,7 +26,7 @@ class StateVectorBase(pybamm.Symbol):
         List of boolean arrays representing slices. Default is None, in which case the
         evaluation_array is computed from y_slices.
 
-    *Extends:* :class:`Array`
+    *Extends:* :class:`pybamm.Symbol`
     """
 
     def __init__(
@@ -212,7 +212,7 @@ class StateVector(StateVectorBase):
         List of boolean arrays representing slices. Default is None, in which case the
         evaluation_array is computed from y_slices.
 
-    *Extends:* :class:`Array`
+    *Extends:* :class:`pybamm.StateVectorBase`
     """
 
     def __init__(
diff --git a/pybamm/expression_tree/variable.py b/pybamm/expression_tree/variable.py
index 3be2f0b61..44055c651 100644
--- a/pybamm/expression_tree/variable.py
+++ b/pybamm/expression_tree/variable.py
@@ -26,20 +26,33 @@ class VariableBase(pybamm.Symbol):
         collector'. For the DFN, the particle concentration would be a Variable with
         domain 'negative particle', secondary domain 'negative electrode' and tertiary
         domain 'current collector'
+    bounds : tuple, optional
+        Physical bounds on the variable
 
     *Extends:* :class:`Symbol`
     """
 
-    def __init__(self, name, domain=None, auxiliary_domains=None):
+    def __init__(self, name, domain=None, auxiliary_domains=None, bounds=None):
         if domain is None:
             domain = []
         if auxiliary_domains is None:
             auxiliary_domains = {}
         super().__init__(name, domain=domain, auxiliary_domains=auxiliary_domains)
+        if bounds is None:
+            bounds = (-np.inf, np.inf)
+        else:
+            if bounds[0] >= bounds[1]:
+                raise ValueError(
+                    "Invalid bounds {}. ".format(bounds)
+                    + "Lower bound should be strictly less than upper bound."
+                )
+        self.bounds = bounds
 
     def new_copy(self):
         """ See :meth:`pybamm.Symbol.new_copy()`. """
-        return self.__class__(self.name, self.domain, self.auxiliary_domains)
+        return self.__class__(
+            self.name, self.domain, self.auxiliary_domains, self.bounds
+        )
 
     def _evaluate_for_shape(self):
         """ See :meth:`pybamm.Symbol.evaluate_for_shape_using_domain()` """
@@ -59,21 +72,24 @@ class Variable(VariableBase):
 
     name : str
         name of the node
-    domain : iterable of str
+        domain : iterable of str, optional
         list of domains that this variable is valid over
-    auxiliary_domains : dict
+    auxiliary_domains : dict, optional
         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'
-
+    bounds : tuple, optional
+        Physical bounds on the variable
     *Extends:* :class:`Symbol`
     """
 
-    def __init__(self, name, domain=None, auxiliary_domains=None):
-        super().__init__(name, domain=domain, auxiliary_domains=auxiliary_domains)
+    def __init__(self, name, domain=None, auxiliary_domains=None, bounds=None):
+        super().__init__(
+            name, domain=domain, auxiliary_domains=auxiliary_domains, bounds=bounds
+        )
 
     def diff(self, variable):
         if variable.id == self.id:
@@ -110,11 +126,13 @@ class VariableDot(VariableBase):
         collector'. For the DFN, the particle concentration would be a Variable with
         domain 'negative particle', secondary domain 'negative electrode' and tertiary
         domain 'current collector'
-
+    bounds : tuple, optional
+        Physical bounds on the variable. Included for compatibility with `VariableBase`,
+        but ignored.
     *Extends:* :class:`Symbol`
     """
 
-    def __init__(self, name, domain=None, auxiliary_domains=None):
+    def __init__(self, name, domain=None, auxiliary_domains=None, bounds=None):
         super().__init__(name, domain=domain, auxiliary_domains=auxiliary_domains)
 
     def get_variable(self):
diff --git a/pybamm/expression_tree/vector.py b/pybamm/expression_tree/vector.py
index 7714ab176..a03f177f4 100644
--- a/pybamm/expression_tree/vector.py
+++ b/pybamm/expression_tree/vector.py
@@ -21,6 +21,8 @@ class Vector(pybamm.Array):
         auxiliary_domains=None,
         entries_string=None,
     ):
+        if isinstance(entries, list):
+            entries = np.array(entries)
         # make sure that entries are a vector (can be a column vector)
         if entries.ndim == 1:
             entries = entries[:, np.newaxis]
diff --git a/pybamm/models/standard_variables.py b/pybamm/models/standard_variables.py
index bf0b99558..87e0c7b03 100644
--- a/pybamm/models/standard_variables.py
+++ b/pybamm/models/standard_variables.py
@@ -2,27 +2,33 @@
 # Standard variables for the models
 #
 import pybamm
+import numpy as np
 
 # Electrolyte concentration
 c_e_n = pybamm.Variable(
     "Negative electrolyte concentration",
     domain="negative electrode",
     auxiliary_domains={"secondary": "current collector"},
+    bounds=(0, np.inf),
 )
 c_e_s = pybamm.Variable(
     "Separator electrolyte concentration",
     domain="separator",
     auxiliary_domains={"secondary": "current collector"},
+    bounds=(0, np.inf),
 )
 c_e_p = pybamm.Variable(
     "Positive electrolyte concentration",
     domain="positive electrode",
     auxiliary_domains={"secondary": "current collector"},
+    bounds=(0, np.inf),
 )
 c_e = pybamm.Concatenation(c_e_n, c_e_s, c_e_p)
 
 c_e_av = pybamm.Variable(
-    "X-averaged electrolyte concentration", domain="current collector"
+    "X-averaged electrolyte concentration",
+    domain="current collector",
+    bounds=(0, np.inf),
 )
 
 # Electrolyte potential
@@ -105,6 +111,7 @@ c_s_n = pybamm.Variable(
         "secondary": "negative electrode",
         "tertiary": "current collector",
     },
+    bounds=(0, 1),
 )
 c_s_p = pybamm.Variable(
     "Positive particle concentration",
@@ -113,32 +120,41 @@ c_s_p = pybamm.Variable(
         "secondary": "positive electrode",
         "tertiary": "current collector",
     },
+    bounds=(0, 1),
 )
 c_s_n_xav = pybamm.Variable(
     "X-averaged negative particle concentration",
     domain="negative particle",
     auxiliary_domains={"secondary": "current collector"},
+    bounds=(0, 1),
 )
 c_s_p_xav = pybamm.Variable(
     "X-averaged positive particle concentration",
     domain="positive particle",
     auxiliary_domains={"secondary": "current collector"},
+    bounds=(0, 1),
 )
 c_s_n_surf = pybamm.Variable(
     "Negative particle surface concentration",
     domain="negative electrode",
     auxiliary_domains={"secondary": "current collector"},
+    bounds=(0, 1),
 )
 c_s_p_surf = pybamm.Variable(
     "Positive particle surface concentration",
     domain="positive electrode",
     auxiliary_domains={"secondary": "current collector"},
+    bounds=(0, 1),
 )
 c_s_n_surf_xav = pybamm.Variable(
-    "X-averaged negative particle surface concentration", domain="current collector"
+    "X-averaged negative particle surface concentration",
+    domain="current collector",
+    bounds=(0, 1),
 )
 c_s_p_surf_xav = pybamm.Variable(
-    "X-averaged positive particle surface concentration", domain="current collector"
+    "X-averaged positive particle surface concentration",
+    domain="current collector",
+    bounds=(0, 1),
 )
 
 
@@ -147,26 +163,31 @@ eps_n = pybamm.Variable(
     "Negative electrode porosity",
     domain="negative electrode",
     auxiliary_domains={"secondary": "current collector"},
+    bounds=(0, 1),
 )
 eps_s = pybamm.Variable(
     "Separator porosity",
     domain="separator",
     auxiliary_domains={"secondary": "current collector"},
+    bounds=(0, 1),
 )
 eps_p = pybamm.Variable(
     "Positive electrode porosity",
     domain="positive electrode",
     auxiliary_domains={"secondary": "current collector"},
+    bounds=(0, 1),
 )
 eps = pybamm.Concatenation(eps_n, eps_s, eps_p)
 
 # Piecewise constant (for asymptotic models)
 eps_n_pc = pybamm.Variable(
-    "X-averaged negative electrode porosity", domain="current collector"
+    "X-averaged negative electrode porosity", domain="current collector", bounds=(0, 1),
+)
+eps_s_pc = pybamm.Variable(
+    "X-averaged separator porosity", domain="current collector", bounds=(0, 1)
 )
-eps_s_pc = pybamm.Variable("X-averaged separator porosity", domain="current collector")
 eps_p_pc = pybamm.Variable(
-    "X-averaged positive electrode porosity", domain="current collector"
+    "X-averaged positive electrode porosity", domain="current collector", bounds=(0, 1),
 )
 
 eps_piecewise_constant = pybamm.Concatenation(
@@ -219,4 +240,3 @@ L_outer = pybamm.Variable(
     domain=["negative electrode"],
     auxiliary_domains={"secondary": "current collector"},
 )
-
diff --git a/pybamm/solvers/algebraic_solver.py b/pybamm/solvers/algebraic_solver.py
index da93e50d9..dd4e0d0b2 100644
--- a/pybamm/solvers/algebraic_solver.py
+++ b/pybamm/solvers/algebraic_solver.py
@@ -19,7 +19,9 @@ class AlgebraicSolver(pybamm.BaseSolver):
     Parameters
     ----------
     method : str, optional
-        The method to use to solve the system (default is "lm")
+        The method to use to solve the system (default is "lm"). If it starts with
+        "lsq", least-squares minimization is used. The method for least-squares can be
+        specified in the form "lsq_methodname"
     tol : float, optional
         The tolerance for the solver (default is 1e-6).
     extra_options : dict, optional
@@ -123,14 +125,68 @@ class AlgebraicSolver(pybamm.BaseSolver):
                 y_alg[:, idx] = y0_alg
             # Otherwise calculate new y0
             else:
-                sol = optimize.root(
-                    root_fun,
-                    y0_alg,
-                    method=self.method,
-                    tol=self.tol,
-                    jac=jac_fn,
-                    options=self.extra_options,
-                )
+                # Methods which use least-squares are specified as either "lsq", which
+                # uses the default method, or with "lsq__methodname"
+                if self.method.startswith("lsq"):
+
+                    if self.method == "lsq":
+                        method = "trf"
+                    else:
+                        method = self.method[5:]
+                    if jac_fn is None:
+                        jac_fn = "2-point"
+                    sol = optimize.least_squares(
+                        root_fun,
+                        y0_alg,
+                        method=method,
+                        ftol=self.tol,
+                        jac=jac_fn,
+                        bounds=model.bounds,
+                        **self.extra_options,
+                    )
+                # Methods which use minimize are specified as either "minimize", which
+                # uses the default method, or with "minimize__methodname"
+                elif self.method.startswith("minimize"):
+                    # Adapt the root function for minimize
+                    def root_norm(y):
+                        return np.sum(root_fun(y) ** 2)
+
+                    if jac_fn is None:
+                        jac_norm = None
+                    else:
+
+                        def jac_norm(y):
+                            return np.sum(2 * root_fun(y) * jac_fn(y), 0)
+
+                    if self.method == "minimize":
+                        method = None
+                    else:
+                        method = self.method[10:]
+                    extra_options = self.extra_options
+                    if np.any(model.bounds[0] != -np.inf) or np.any(
+                        model.bounds[1] != np.inf
+                    ):
+                        bounds = [
+                            (lb, ub) for lb, ub in zip(model.bounds[0], model.bounds[1])
+                        ]
+                        extra_options["bounds"] = bounds
+                    sol = optimize.minimize(
+                        root_norm,
+                        y0_alg,
+                        method=method,
+                        tol=self.tol,
+                        jac=jac_norm,
+                        **extra_options,
+                    )
+                else:
+                    sol = optimize.root(
+                        root_fun,
+                        y0_alg,
+                        method=self.method,
+                        tol=self.tol,
+                        jac=jac_fn,
+                        options=self.extra_options,
+                    )
 
                 if sol.success and np.all(abs(sol.fun) < self.tol):
                     # update initial guess for the next iteration
@@ -143,12 +199,10 @@ class AlgebraicSolver(pybamm.BaseSolver):
                     )
                 else:
                     raise pybamm.SolverError(
-                        """
-                        Could not find acceptable solution: solver terminated
-                        successfully, but maximum solution error ({})
-                        above tolerance ({})
-                        """.format(
-                            np.max(sol.fun), self.tol
+                        "Could not find acceptable solution: solver terminated "
+                        "successfully, but maximum solution error "
+                        "({}) above tolerance ({})".format(
+                            np.max(abs(sol.fun)), self.tol
                         )
                     )
 
diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py
index 159118bf7..321709781 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -51,7 +51,7 @@ class BaseSolver(object):
                 "max_steps has been deprecated, and should be set using the "
                 "solver-specific extra-options dictionaries instead"
             )
-        self.models_set_up = set()
+        self.models_set_up = {}
 
         # Defaults, can be overwritten by specific solver
         self.name = "Base solver"
@@ -114,7 +114,7 @@ class BaseSolver(object):
         "Returns a copy of the solver"
         new_solver = copy.copy(self)
         # clear models_set_up
-        new_solver.models_set_up = set()
+        new_solver.models_set_up = {}
         return new_solver
 
     def set_up(self, model, inputs=None):
@@ -504,9 +504,23 @@ class BaseSolver(object):
         if model not in self.models_set_up:
             self.set_up(model, ext_and_inputs)
             set_up_time = timer.time()
-            self.models_set_up.add(model)
+            self.models_set_up.update(
+                {model: {"initial conditions": model.concatenated_initial_conditions}}
+            )
         else:
-            set_up_time = 0
+            ics_set_up = self.models_set_up[model]["initial conditions"]
+            # Check that initial conditions have not been updated
+            if ics_set_up.id == model.concatenated_initial_conditions.id:
+                set_up_time = 0
+            else:
+                # If the new initial conditions are different, set up again
+                # Doing the whole setup again might be slow, but no need to prematurely
+                # optimize this
+                self.set_up(model, ext_and_inputs)
+                self.models_set_up[model][
+                    "initial conditions"
+                ] = model.concatenated_initial_conditions
+                set_up_time = timer.time()
 
         # (Re-)calculate consistent initial conditions
         self._set_initial_conditions(model, ext_and_inputs, update_rhs=True)
diff --git a/pybamm/solvers/casadi_algebraic_solver.py b/pybamm/solvers/casadi_algebraic_solver.py
index 229cc7f6d..e286936f2 100644
--- a/pybamm/solvers/casadi_algebraic_solver.py
+++ b/pybamm/solvers/casadi_algebraic_solver.py
@@ -68,6 +68,7 @@ class CasadiAlgebraicSolver(pybamm.BaseSolver):
         # equations will be equal to the initial condition provided. This allows this
         # solver to be used for initialising the DAE solvers
         if model.rhs == {}:
+            len_rhs = 0
             y0_diff = casadi.DM()
             y0_alg = y0
         else:
@@ -86,12 +87,25 @@ class CasadiAlgebraicSolver(pybamm.BaseSolver):
         t_p_sym = casadi.vertcat(t_sym, p_sym)
         alg = model.casadi_algebraic(t_sym, y_sym, p_sym)
 
+        # Set constraints vector in the casadi format
+        # Constrain the unknowns. 0 (default): no constraint on ui, 1: ui >= 0.0,
+        # -1: ui <= 0.0, 2: ui > 0.0, -2: ui < 0.0.
+        constraints = np.zeros_like(model.bounds[0], dtype=int)
+        # If the lower bound is positive then the variable must always be positive
+        constraints[model.bounds[0] >= 0] = 1
+        # If the upper bound is negative then the variable must always be negative
+        constraints[model.bounds[1] <= 0] = -1
+
         # Set up rootfinder
         roots = casadi.rootfinder(
             "roots",
             "newton",
             dict(x=y_alg_sym, p=t_p_sym, g=alg),
-            {**self.extra_options, "abstol": self.tol},
+            {
+                **self.extra_options,
+                "abstol": self.tol,
+                "constraints": list(constraints[len_rhs:]),
+            },
         )
         for idx, t in enumerate(t_eval):
             # Evaluate algebraic with new t and previous y0, if it's already close
@@ -146,7 +160,7 @@ class CasadiAlgebraicSolver(pybamm.BaseSolver):
                         successfully, but maximum solution error ({})
                         above tolerance ({})
                         """.format(
-                            casadi.mmax(fun), self.tol
+                            casadi.mmax(casadi.fabs(fun)), self.tol
                         )
                     )
 
diff --git a/pybamm/spatial_methods/finite_volume.py b/pybamm/spatial_methods/finite_volume.py
index c2e9683cb..8b00a6bfc 100644
--- a/pybamm/spatial_methods/finite_volume.py
+++ b/pybamm/spatial_methods/finite_volume.py
@@ -641,7 +641,7 @@ class FiniteVolume(pybamm.SpatialMethod):
         n = submesh.npts
         second_dim_repeats = self._get_auxiliary_domain_repeats(symbol.domains)
 
-        bcs_vector = pybamm.Vector(np.array([]))  # starts empty
+        bcs_vector = pybamm.Vector([])  # starts empty
 
         lbc_value, lbc_type = bcs["left"]
         rbc_value, rbc_type = bcs["right"]

Adding constraints to algebraic solver
It seems the algebraic solver sometimes fails when a state goes outside an obvious bound (e.g. concentration going negative) which makes the algebraic equations return nan. Opening this issue to brainstorm some solutions, so suggestions welcome:

1. Simply replace `nan` with a very large number. Will try this first and hope it works as this would be the simplest solution. The other two solutions require manually specifying some bounds in the models (e.g. using the existing event functionality)
2. Impose hard constraints, but neither `scipy.optimize.root` nor casadi have a very clean way of providing constraints
3. Impose soft constraints: add a function of `c` with singularities as c-> 0 and c -> 1 to the function being optimized

N.B. the use case here is for fuel cell models, solving directly for the steady state. 
diff --git a/tests/unit/test_discretisations/test_discretisation.py b/tests/unit/test_discretisations/test_discretisation.py
index f827a5a13..7ceff8746 100644
--- a/tests/unit/test_discretisations/test_discretisation.py
+++ b/tests/unit/test_discretisations/test_discretisation.py
@@ -131,6 +131,7 @@ class TestDiscretise(unittest.TestCase):
         self.assertEqual(disc.y_slices[a.id][0], slice(0, 10, None))
 
         self.assertEqual(model.y_slices[a][0], slice(0, 10, None))
+        self.assertEqual(model.bounds, disc.bounds)
 
         b_test = np.ones((10, 1))
         np.testing.assert_array_equal(
@@ -265,7 +266,7 @@ class TestDiscretise(unittest.TestCase):
         np.testing.assert_array_equal(y[disc.y_slices[c.id][0]], c_true)
 
         # Several variables
-        d = pybamm.Variable("d", domain=whole_cell)
+        d = pybamm.Variable("d", domain=whole_cell, bounds=(0, 1))
         jn = pybamm.Variable("jn", domain=["negative electrode"])
         variables = [c, d, jn]
         disc.set_variable_slices(variables)
@@ -274,6 +275,12 @@ class TestDiscretise(unittest.TestCase):
             disc.y_slices,
             {c.id: [slice(0, 100)], d.id: [slice(100, 200)], jn.id: [slice(200, 240)]},
         )
+        np.testing.assert_array_equal(
+            disc.bounds[0], [-np.inf] * 100 + [0] * 100 + [-np.inf] * 40
+        )
+        np.testing.assert_array_equal(
+            disc.bounds[1], [np.inf] * 100 + [1] * 100 + [np.inf] * 40
+        )
         d_true = 4 * combined_submesh.nodes
         jn_true = mesh["negative electrode"].nodes ** 3
         y = np.concatenate([c_true, d_true, jn_true])
@@ -297,6 +304,12 @@ class TestDiscretise(unittest.TestCase):
                 jp.id: [slice(265, 300)],
             },
         )
+        np.testing.assert_array_equal(
+            disc.bounds[0], [-np.inf] * 100 + [0] * 100 + [-np.inf] * 100
+        )
+        np.testing.assert_array_equal(
+            disc.bounds[1], [np.inf] * 100 + [1] * 100 + [np.inf] * 100
+        )
         d_true = 4 * combined_submesh.nodes
         jn_true = mesh["negative electrode"].nodes ** 3
         y = np.concatenate([c_true, d_true, jn_true])
@@ -338,12 +351,12 @@ class TestDiscretise(unittest.TestCase):
         self.assertIsInstance(scal_disc, pybamm.Scalar)
         self.assertEqual(scal_disc.value, scal.value)
         # vector
-        vec = pybamm.Vector(np.array([1, 2, 3, 4]))
+        vec = pybamm.Vector([1, 2, 3, 4])
         vec_disc = disc.process_symbol(vec)
         self.assertIsInstance(vec_disc, pybamm.Vector)
         np.testing.assert_array_equal(vec_disc.entries, vec.entries)
         # matrix
-        mat = pybamm.Matrix(np.array([[1, 2, 3, 4], [5, 6, 7, 8]]))
+        mat = pybamm.Matrix([[1, 2, 3, 4], [5, 6, 7, 8]])
         mat_disc = disc.process_symbol(mat)
         self.assertIsInstance(mat_disc, pybamm.Matrix)
         np.testing.assert_array_equal(mat_disc.entries, mat.entries)
@@ -839,39 +852,39 @@ class TestDiscretise(unittest.TestCase):
         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"])
-    #     cs = pybamm.Variable("c", domain=["separator"])
-    #     cp = pybamm.Variable("c", domain=["positive electrode"])
-    #     c = pybamm.Concatenation(cn, cs, cp)
-    #     N = pybamm.grad(c)
-    #     model = pybamm.BaseModel()
-    #     model.rhs = {c: pybamm.div(N)}
-    #     model.initial_conditions = {c: pybamm.Scalar(3)}
-
-    #     model.boundary_conditions = {
-    #         c: {"left": (0, "Neumann"), "right": (0, "Neumann")}
-    #     }
-    #     model.check_well_posedness()
-
-    #     # create discretisation
-    #     disc = get_discretisation_for_testing()
-    #     mesh = disc.mesh
-
-    #     combined_submesh = mesh.combine_submeshes(
-    #         "negative electrode", "separator", "positive electrode"
-    #     )
-
-    #     disc.process_model(model)
-    #     y0 = model.concatenated_initial_conditions.evaluate()
-    #     np.testing.assert_array_equal(
-    #         y0, 3 * np.ones_like(combined_submesh.nodes[:, np.newaxis])
-    #     )
-
-    #     # grad and div are identity operators here
-    #     np.testing.assert_array_equal(y0, model.concatenated_rhs.evaluate(None, y0))
-    #     model.check_well_posedness()
+    def test_process_model_concatenation(self):
+        # concatenation of variables as the key
+        cn = pybamm.Variable("c", domain=["negative electrode"])
+        cs = pybamm.Variable("c", domain=["separator"])
+        cp = pybamm.Variable("c", domain=["positive electrode"])
+        c = pybamm.Concatenation(cn, cs, cp)
+        N = pybamm.grad(c)
+        model = pybamm.BaseModel()
+        model.rhs = {c: pybamm.div(N)}
+        model.initial_conditions = {c: pybamm.Scalar(3)}
+
+        model.boundary_conditions = {
+            c: {"left": (0, "Neumann"), "right": (0, "Neumann")}
+        }
+        model.check_well_posedness()
+
+        # create discretisation
+        disc = get_discretisation_for_testing()
+        mesh = disc.mesh
+
+        combined_submesh = mesh.combine_submeshes(
+            "negative electrode", "separator", "positive electrode"
+        )
+
+        disc.process_model(model)
+        y0 = model.concatenated_initial_conditions.evaluate()
+        np.testing.assert_array_equal(
+            y0, 3 * np.ones_like(combined_submesh.nodes[:, np.newaxis])
+        )
+
+        # grad and div are identity operators here
+        np.testing.assert_array_equal(y0, model.concatenated_rhs.evaluate(None, y0))
+        model.check_well_posedness()
 
     def test_process_model_not_inplace(self):
         # concatenation of variables as the key
@@ -1060,16 +1073,19 @@ class TestDiscretise(unittest.TestCase):
             "a",
             domain=["negative electrode"],
             auxiliary_domains={"secondary": "current collector"},
+            bounds=(-5, -2),
         )
         b = pybamm.Variable(
             "b",
             domain=["separator"],
             auxiliary_domains={"secondary": "current collector"},
+            bounds=(6, 10),
         )
         c = pybamm.Variable(
             "c",
             domain=["positive electrode"],
             auxiliary_domains={"secondary": "current collector"},
+            bounds=(0, 1),
         )
         conc = pybamm.Concatenation(a, b, c)
         disc.set_variable_slices([conc])
@@ -1082,6 +1098,12 @@ class TestDiscretise(unittest.TestCase):
         self.assertEqual(
             disc.y_slices[c.id], [slice(65, 100), slice(165, 200), slice(265, 300)]
         )
+        np.testing.assert_array_equal(
+            disc.bounds[0], ([-5] * 40 + [6] * 25 + [0] * 35) * 3
+        )
+        np.testing.assert_array_equal(
+            disc.bounds[1], ([-2] * 40 + [10] * 25 + [1] * 35) * 3
+        )
         expr = disc.process_symbol(conc)
         self.assertIsInstance(expr, pybamm.DomainConcatenation)
 
diff --git a/tests/unit/test_expression_tree/test_array.py b/tests/unit/test_expression_tree/test_array.py
index 550abd84c..e2adab03b 100644
--- a/tests/unit/test_expression_tree/test_array.py
+++ b/tests/unit/test_expression_tree/test_array.py
@@ -12,6 +12,12 @@ class TestArray(unittest.TestCase):
         arr = pybamm.Array(np.array([1, 2, 3]))
         self.assertEqual(arr.name, "Array of shape (3, 1)")
 
+    def test_list_entries(self):
+        vect = pybamm.Array([1, 2, 3])
+        np.testing.assert_array_equal(vect.entries, np.array([[1], [2], [3]]))
+        vect = pybamm.Array([[1], [2], [3]])
+        np.testing.assert_array_equal(vect.entries, np.array([[1], [2], [3]]))
+
     def test_linspace(self):
         x = np.linspace(0, 1, 100)[:, np.newaxis]
         y = pybamm.linspace(0, 1, 100)
diff --git a/tests/unit/test_expression_tree/test_binary_operators.py b/tests/unit/test_expression_tree/test_binary_operators.py
index 02a4b5cea..71d1689c3 100644
--- a/tests/unit/test_expression_tree/test_binary_operators.py
+++ b/tests/unit/test_expression_tree/test_binary_operators.py
@@ -336,7 +336,7 @@ class TestIsZero(unittest.TestCase):
 
         a = pybamm.Matrix(np.zeros((10, 10)))
         b = pybamm.Matrix(np.ones((10, 10)))
-        c = pybamm.Matrix(np.array([1, 0, 0]))
+        c = pybamm.Matrix([1, 0, 0])
         self.assertTrue(pybamm.is_matrix_zero(a))
         self.assertFalse(pybamm.is_matrix_zero(b))
         self.assertFalse(pybamm.is_matrix_zero(c))
diff --git a/tests/unit/test_expression_tree/test_concatenations.py b/tests/unit/test_expression_tree/test_concatenations.py
index 3e2e42b2b..e87f2aaa2 100644
--- a/tests/unit/test_expression_tree/test_concatenations.py
+++ b/tests/unit/test_expression_tree/test_concatenations.py
@@ -18,9 +18,9 @@ class TestConcatenations(unittest.TestCase):
         self.assertEqual(conc.children[0].name, "a")
         self.assertEqual(conc.children[1].name, "b")
         self.assertEqual(conc.children[2].name, "c")
-        d = pybamm.Vector(np.array([2]))
-        e = pybamm.Vector(np.array([1]))
-        f = pybamm.Vector(np.array([3]))
+        d = pybamm.Vector([2])
+        e = pybamm.Vector([1])
+        f = pybamm.Vector([3])
         conc2 = pybamm.Concatenation(d, e, f)
         with self.assertRaises(TypeError):
             conc2.evaluate()
diff --git a/tests/unit/test_expression_tree/test_matrix.py b/tests/unit/test_expression_tree/test_matrix.py
index 0d2300ed5..f4eb8988d 100644
--- a/tests/unit/test_expression_tree/test_matrix.py
+++ b/tests/unit/test_expression_tree/test_matrix.py
@@ -19,6 +19,12 @@ class TestMatrix(unittest.TestCase):
         self.assertEqual(self.mat.shape, (3, 3))
         self.assertEqual(self.mat.size, 9)
 
+    def test_list_entry(self):
+        mat = pybamm.Matrix([[1, 2, 0], [0, 1, 0], [0, 0, 1]])
+        np.testing.assert_array_equal(
+            mat.entries, np.array([[1, 2, 0], [0, 1, 0], [0, 0, 1]])
+        )
+
     def test_matrix_evaluate(self):
         np.testing.assert_array_equal(self.mat.evaluate(), self.A)
 
diff --git a/tests/unit/test_expression_tree/test_operations/test_copy.py b/tests/unit/test_expression_tree/test_operations/test_copy.py
index ef6d3a063..564eee5e6 100644
--- a/tests/unit/test_expression_tree/test_operations/test_copy.py
+++ b/tests/unit/test_expression_tree/test_operations/test_copy.py
@@ -14,7 +14,7 @@ class TestCopy(unittest.TestCase):
         v_n = pybamm.Variable("v", "negative electrode")
         x_n = pybamm.standard_spatial_vars.x_n
         v_s = pybamm.Variable("v", "separator")
-        vec = pybamm.Vector(np.array([1, 2, 3, 4, 5]))
+        vec = pybamm.Vector([1, 2, 3, 4, 5])
         mesh = get_mesh_for_testing()
 
         for symbol in [
diff --git a/tests/unit/test_expression_tree/test_operations/test_evaluate.py b/tests/unit/test_expression_tree/test_operations/test_evaluate.py
index 8ff7eaaf4..3f08f1137 100644
--- a/tests/unit/test_expression_tree/test_operations/test_evaluate.py
+++ b/tests/unit/test_expression_tree/test_operations/test_evaluate.py
@@ -106,7 +106,7 @@ class TestEvaluate(unittest.TestCase):
         # test matrix
         constant_symbols = OrderedDict()
         variable_symbols = OrderedDict()
-        A = pybamm.Matrix(np.array([[1, 2], [3, 4]]))
+        A = pybamm.Matrix([[1, 2], [3, 4]])
         pybamm.find_symbols(A, constant_symbols, variable_symbols)
         self.assertEqual(len(variable_symbols), 0)
         self.assertEqual(list(constant_symbols.keys())[0], A.id)
@@ -335,7 +335,7 @@ class TestEvaluate(unittest.TestCase):
             self.assertEqual(result, expr.evaluate(t=t, y=y))
 
         # test something with a matrix multiplication
-        A = pybamm.Matrix(np.array([[1, 2], [3, 4]]))
+        A = pybamm.Matrix([[1, 2], [3, 4]])
         expr = A @ pybamm.StateVector(slice(0, 2))
         evaluator = pybamm.EvaluatorPython(expr)
         for t, y in zip(t_tests, y_tests):
@@ -343,7 +343,7 @@ class TestEvaluate(unittest.TestCase):
             np.testing.assert_allclose(result, expr.evaluate(t=t, y=y))
 
         # test something with a heaviside
-        a = pybamm.Vector(np.array([1, 2]))
+        a = pybamm.Vector([1, 2])
         expr = a <= pybamm.StateVector(slice(0, 2))
         evaluator = pybamm.EvaluatorPython(expr)
         for t, y in zip(t_tests, y_tests):
@@ -357,7 +357,7 @@ class TestEvaluate(unittest.TestCase):
             np.testing.assert_allclose(result, expr.evaluate(t=t, y=y))
 
         # test something with a minimum or maximum
-        a = pybamm.Vector(np.array([1, 2]))
+        a = pybamm.Vector([1, 2])
         expr = pybamm.minimum(a, pybamm.StateVector(slice(0, 2)))
         evaluator = pybamm.EvaluatorPython(expr)
         for t, y in zip(t_tests, y_tests):
@@ -378,7 +378,7 @@ class TestEvaluate(unittest.TestCase):
             self.assertEqual(result, expr.evaluate(t=t, y=y))
 
         # test something with a sparse matrix multiplication
-        A = pybamm.Matrix(np.array([[1, 2], [3, 4]]))
+        A = pybamm.Matrix([[1, 2], [3, 4]])
         B = pybamm.Matrix(scipy.sparse.csr_matrix(np.array([[1, 0], [0, 4]])))
         C = pybamm.Matrix(scipy.sparse.coo_matrix(np.array([[1, 0], [0, 4]])))
         expr = A @ B @ C @ pybamm.StateVector(slice(0, 2))
@@ -388,8 +388,8 @@ class TestEvaluate(unittest.TestCase):
             np.testing.assert_allclose(result, expr.evaluate(t=t, y=y))
 
         # test numpy concatenation
-        a = pybamm.Vector(np.array([[1], [2]]))
-        b = pybamm.Vector(np.array([[3]]))
+        a = pybamm.Vector([[1], [2]])
+        b = pybamm.Vector([[3]])
         expr = pybamm.NumpyConcatenation(a, b)
         evaluator = pybamm.EvaluatorPython(expr)
         for t, y in zip(t_tests, y_tests):
diff --git a/tests/unit/test_expression_tree/test_operations/test_jac_2D.py b/tests/unit/test_expression_tree/test_operations/test_jac_2D.py
index 3bd80db91..5fea85ee7 100644
--- a/tests/unit/test_expression_tree/test_operations/test_jac_2D.py
+++ b/tests/unit/test_expression_tree/test_operations/test_jac_2D.py
@@ -222,6 +222,10 @@ class TestJacobian(unittest.TestCase):
         a = 2 * pybamm.PrimaryBroadcast(curr_coll_vector, a_dom)
         b = pybamm.PrimaryBroadcast(curr_coll_vector, b_dom)
         c = 3 * pybamm.PrimaryBroadcast(curr_coll_vector, c_dom)
+        # Add bounds for compatibility with the discretisation
+        a.bounds = (-np.inf, np.inf)
+        b.bounds = (-np.inf, np.inf)
+        c.bounds = (-np.inf, np.inf)
 
         conc = pybamm.Concatenation(a, b, c)
         disc.set_variable_slices([conc])
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 b74259d89..3aac02b6f 100644
--- a/tests/unit/test_expression_tree/test_operations/test_simplify.py
+++ b/tests/unit/test_expression_tree/test_operations/test_simplify.py
@@ -336,7 +336,7 @@ class TestSimplify(unittest.TestCase):
         b = pybamm.Matrix(np.ones((2, 2)))
 
         # matrix multiplication
-        A = pybamm.Matrix(np.array([[1, 0], [0, 1]]))
+        A = pybamm.Matrix([[1, 0], [0, 1]])
         self.assertIsInstance((a @ A).simplify(), pybamm.Matrix)
         np.testing.assert_array_equal(
             (a @ A).simplify().evaluate().toarray(), np.zeros((2, 2))
@@ -347,8 +347,8 @@ class TestSimplify(unittest.TestCase):
         )
 
         # matrix * matrix
-        m1 = pybamm.Matrix(np.array([[2, 0], [0, 2]]))
-        m2 = pybamm.Matrix(np.array([[3, 0], [0, 3]]))
+        m1 = pybamm.Matrix([[2, 0], [0, 2]])
+        m2 = pybamm.Matrix([[3, 0], [0, 3]])
         v = pybamm.StateVector(slice(0, 2))
         v2 = pybamm.StateVector(slice(2, 4))
 
@@ -412,8 +412,8 @@ class TestSimplify(unittest.TestCase):
             )
 
         # matrix * vector
-        m1 = pybamm.Matrix(np.array([[2, 0], [0, 2]]))
-        v1 = pybamm.Vector(np.array([1, 1]))
+        m1 = pybamm.Matrix([[2, 0], [0, 2]])
+        v1 = pybamm.Vector([1, 1])
 
         for expr in [(m1 @ v1).simplify()]:
             self.assertIsInstance(expr, pybamm.Vector)
diff --git a/tests/unit/test_expression_tree/test_variable.py b/tests/unit/test_expression_tree/test_variable.py
index 8e4a8be11..b023f1e43 100644
--- a/tests/unit/test_expression_tree/test_variable.py
+++ b/tests/unit/test_expression_tree/test_variable.py
@@ -33,6 +33,18 @@ class TestVariable(unittest.TestCase):
         self.assertNotEqual(a1.id, a3.id)
         self.assertNotEqual(a1.id, a4.id)
 
+    def test_variable_bounds(self):
+        var = pybamm.Variable("var")
+        self.assertEqual(var.bounds, (-np.inf, np.inf))
+
+        var = pybamm.Variable("var", bounds=(0, 1))
+        self.assertEqual(var.bounds, (0, 1))
+
+        with self.assertRaisesRegex(ValueError, "Invalid bounds"):
+            pybamm.Variable("var", bounds=(1, 0))
+        with self.assertRaisesRegex(ValueError, "Invalid bounds"):
+            pybamm.Variable("var", bounds=(1, 1))
+
 
 class TestVariableDot(unittest.TestCase):
     def test_variable_init(self):
diff --git a/tests/unit/test_expression_tree/test_vector.py b/tests/unit/test_expression_tree/test_vector.py
index 924876dd8..52c414b78 100644
--- a/tests/unit/test_expression_tree/test_vector.py
+++ b/tests/unit/test_expression_tree/test_vector.py
@@ -21,6 +21,12 @@ class TestVector(unittest.TestCase):
         vect1d = pybamm.Vector(np.array([1, 2, 3]))
         np.testing.assert_array_equal(self.vect.entries, vect1d.entries)
 
+    def test_list_entries(self):
+        vect = pybamm.Vector([1, 2, 3])
+        np.testing.assert_array_equal(vect.entries, np.array([[1], [2], [3]]))
+        vect = pybamm.Vector([[1], [2], [3]])
+        np.testing.assert_array_equal(vect.entries, np.array([[1], [2], [3]]))
+
     def test_vector_evaluate(self):
         np.testing.assert_array_equal(self.vect.evaluate(), self.x)
 
diff --git a/tests/unit/test_solvers/test_algebraic_solver.py b/tests/unit/test_solvers/test_algebraic_solver.py
index 45c858327..4e7d660b5 100644
--- a/tests/unit/test_solvers/test_algebraic_solver.py
+++ b/tests/unit/test_solvers/test_algebraic_solver.py
@@ -135,6 +135,7 @@ class TestAlgebraicSolver(unittest.TestCase):
 
         # Test without jacobian
         model.use_jacobian = False
+        solver.models_set_up = set()
         solution_no_jac = solver.solve(model)
         np.testing.assert_array_equal(
             model.variables["var1"].evaluate(t=None, y=solution_no_jac.y), sol[:100]
@@ -143,6 +144,112 @@ class TestAlgebraicSolver(unittest.TestCase):
             model.variables["var2"].evaluate(t=None, y=solution_no_jac.y), sol[100:]
         )
 
+    def test_model_solver_least_squares(self):
+        # Create model
+        model = pybamm.BaseModel()
+        whole_cell = ["negative electrode", "separator", "positive electrode"]
+        var1 = pybamm.Variable("var1", domain=whole_cell)
+        var2 = pybamm.Variable("var2", domain=whole_cell)
+        model.algebraic = {var1: var1 - 3, var2: 2 * var1 - var2}
+        model.initial_conditions = {var1: pybamm.Scalar(1), var2: pybamm.Scalar(4)}
+        model.variables = {"var1": var1, "var2": var2}
+        disc = get_discretisation_for_testing()
+        disc.process_model(model)
+
+        sol = np.concatenate((np.ones(100) * 3, np.ones(100) * 6))[:, np.newaxis]
+
+        # Solve
+        solver = pybamm.AlgebraicSolver("lsq")
+        solution = solver.solve(model)
+        np.testing.assert_array_almost_equal(
+            model.variables["var1"].evaluate(t=None, y=solution.y), sol[:100]
+        )
+        np.testing.assert_array_almost_equal(
+            model.variables["var2"].evaluate(t=None, y=solution.y), sol[100:]
+        )
+
+        # Test without jacobian and with a different method
+        model.use_jacobian = False
+        solver = pybamm.AlgebraicSolver("lsq__trf")
+        solution_no_jac = solver.solve(model)
+        np.testing.assert_array_almost_equal(
+            model.variables["var1"].evaluate(t=None, y=solution_no_jac.y), sol[:100]
+        )
+        np.testing.assert_array_almost_equal(
+            model.variables["var2"].evaluate(t=None, y=solution_no_jac.y), sol[100:]
+        )
+
+    def test_model_solver_minimize(self):
+        # Create model
+        model = pybamm.BaseModel()
+        whole_cell = ["negative electrode", "separator", "positive electrode"]
+        var1 = pybamm.Variable("var1", domain=whole_cell)
+        var2 = pybamm.Variable("var2", domain=whole_cell)
+        model.algebraic = {var1: var1 - 3, var2: 2 * var1 - var2}
+        model.initial_conditions = {var1: pybamm.Scalar(1), var2: pybamm.Scalar(4)}
+        model.variables = {"var1": var1, "var2": var2}
+        disc = get_discretisation_for_testing()
+        disc.process_model(model)
+
+        sol = np.concatenate((np.ones(100) * 3, np.ones(100) * 6))[:, np.newaxis]
+
+        # Solve
+        solver = pybamm.AlgebraicSolver("minimize", tol=1e-8)
+        solution = solver.solve(model)
+        np.testing.assert_array_almost_equal(
+            model.variables["var1"].evaluate(t=None, y=solution.y), sol[:100]
+        )
+        np.testing.assert_array_almost_equal(
+            model.variables["var2"].evaluate(t=None, y=solution.y), sol[100:]
+        )
+
+        # Test without jacobian and with a different method
+        model.use_jacobian = False
+        solver = pybamm.AlgebraicSolver("minimize__BFGS")
+        solution_no_jac = solver.solve(model)
+        np.testing.assert_array_almost_equal(
+            model.variables["var1"].evaluate(t=None, y=solution_no_jac.y), sol[:100]
+        )
+        np.testing.assert_array_almost_equal(
+            model.variables["var2"].evaluate(t=None, y=solution_no_jac.y), sol[100:]
+        )
+
+    def test_model_solver_least_squares_with_bounds(self):
+        # Note: we need a better test case to test this functionality properly
+        # Create model
+        model = pybamm.BaseModel()
+        var1 = pybamm.Variable("var1", bounds=(0, 10))
+        model.algebraic = {var1: pybamm.sin(var1) + 1}
+        model.initial_conditions = {var1: pybamm.Scalar(3)}
+        model.variables = {"var1": var1}
+
+        # Solve
+        solver = pybamm.AlgebraicSolver("lsq", tol=1e-5)
+        solution = solver.solve(model)
+        np.testing.assert_array_almost_equal(
+            model.variables["var1"].evaluate(t=None, y=solution.y),
+            3 * np.pi / 2,
+            decimal=2,
+        )
+
+    def test_model_solver_minimize_with_bounds(self):
+        # Note: we need a better test case to test this functionality properly
+        # Create model
+        model = pybamm.BaseModel()
+        var1 = pybamm.Variable("var1", bounds=(0, 10))
+        model.algebraic = {var1: pybamm.sin(var1) + 1}
+        model.initial_conditions = {var1: pybamm.Scalar(3)}
+        model.variables = {"var1": var1}
+
+        # Solve
+        solver = pybamm.AlgebraicSolver("minimize", tol=1e-16)
+        solution = solver.solve(model)
+        np.testing.assert_array_almost_equal(
+            model.variables["var1"].evaluate(t=None, y=solution.y),
+            3 * np.pi / 2,
+            decimal=4,
+        )
+
     def test_model_solver_with_time(self):
         # Create model
         model = pybamm.BaseModel()
diff --git a/tests/unit/test_solvers/test_base_solver.py b/tests/unit/test_solvers/test_base_solver.py
index 3633cf6bf..35e6cf25d 100644
--- a/tests/unit/test_solvers/test_base_solver.py
+++ b/tests/unit/test_solvers/test_base_solver.py
@@ -117,6 +117,7 @@ class TestBaseSolver(unittest.TestCase):
                     "alg", [t, y, p], [self.algebraic_eval(t, y, p)]
                 )
                 self.convert_to_format = "casadi"
+                self.bounds = (np.array([-np.inf]), np.array([np.inf]))
 
             def rhs_eval(self, t, y, inputs):
                 return np.array([])
@@ -151,6 +152,7 @@ class TestBaseSolver(unittest.TestCase):
                     "alg", [t, y, p], [self.algebraic_eval(t, y, p)]
                 )
                 self.convert_to_format = "casadi"
+                self.bounds = (-np.inf * np.ones(4), np.inf * np.ones(4))
 
             def rhs_eval(self, t, y, inputs):
                 return y[0:1]
@@ -197,6 +199,7 @@ class TestBaseSolver(unittest.TestCase):
                     "alg", [t, y, p], [self.algebraic_eval(t, y, p)]
                 )
                 self.convert_to_format = "casadi"
+                self.bounds = (np.array([-np.inf]), np.array([np.inf]))
 
             def rhs_eval(self, t, y, inputs):
                 return np.array([])
diff --git a/tests/unit/test_solvers/test_casadi_algebraic_solver.py b/tests/unit/test_solvers/test_casadi_algebraic_solver.py
index 208860d04..1d8422318 100644
--- a/tests/unit/test_solvers/test_casadi_algebraic_solver.py
+++ b/tests/unit/test_solvers/test_casadi_algebraic_solver.py
@@ -58,6 +58,7 @@ class TestCasadiAlgebraicSolver(unittest.TestCase):
             p = casadi.MX.sym("p")
             rhs = {}
             casadi_algebraic = casadi.Function("alg", [t, y, p], [y ** 2 + 1])
+            bounds = (np.array([-np.inf]), np.array([np.inf]))
 
             def algebraic_eval(self, t, y, inputs):
                 # algebraic equation has no real root
@@ -104,6 +105,22 @@ class TestCasadiAlgebraicSolver(unittest.TestCase):
             sol[1, :],
         )
 
+    def test_model_solver_with_bounds(self):
+        # Note: we need a better test case to test this functionality properly
+        # Create model
+        model = pybamm.BaseModel()
+        var1 = pybamm.Variable("var1", bounds=(0, 10))
+        model.algebraic = {var1: pybamm.sin(var1) + 1}
+        model.initial_conditions = {var1: pybamm.Scalar(3)}
+        model.variables = {"var1": var1}
+
+        # Solve
+        solver = pybamm.CasadiAlgebraicSolver(tol=1e-12)
+        solution = solver.solve(model)
+        np.testing.assert_array_almost_equal(
+            model.variables["var1"].evaluate(t=None, y=solution.y), 3 * np.pi / 2
+        )
+
     def test_solve_with_input(self):
         # Simple system: a single algebraic equation
         var = pybamm.Variable("var")
diff --git a/tests/unit/test_solvers/test_scikits_solvers.py b/tests/unit/test_solvers/test_scikits_solvers.py
index c984d7f40..6ab8fe2dc 100644
--- a/tests/unit/test_solvers/test_scikits_solvers.py
+++ b/tests/unit/test_solvers/test_scikits_solvers.py
@@ -42,7 +42,7 @@ class TestScikitsSolvers(unittest.TestCase):
 
         # Create custom model so that custom ics
         class Model:
-            mass_matrix = pybamm.Matrix(np.array([[1.0, 0.0], [0.0, 0.0]]))
+            mass_matrix = pybamm.Matrix([[1.0, 0.0], [0.0, 0.0]])
             y0 = np.array([0.0, 1.0])
             terminate_events_eval = []
             timescale_eval = 1
@@ -92,7 +92,7 @@ class TestScikitsSolvers(unittest.TestCase):
 
         # Create custom model so that custom mass matrix can be used
         class Model:
-            mass_matrix = pybamm.Matrix(np.array([[4.0, 0.0], [0.0, 0.0]]))
+            mass_matrix = pybamm.Matrix([[4.0, 0.0], [0.0, 0.0]])
             y0 = np.array([0.0, 0.0])
             terminate_events_eval = []
             timescale_eval = 1
diff --git a/tests/unit/test_solvers/test_scipy_solver.py b/tests/unit/test_solvers/test_scipy_solver.py
index 1d6cfd43a..1d0aae34b 100644
--- a/tests/unit/test_solvers/test_scipy_solver.py
+++ b/tests/unit/test_solvers/test_scipy_solver.py
@@ -314,6 +314,30 @@ class TestScipySolver(unittest.TestCase):
             solution.y[0], 1 * np.exp(-0.1 * solution.t), decimal=5
         )
 
+    def test_model_solver_manually_update_initial_conditions(self):
+        # Create model
+        model = pybamm.BaseModel()
+        var1 = pybamm.Variable("var1")
+        model.rhs = {var1: -var1}
+        model.initial_conditions = {var1: 1}
+
+        # Solve
+        solver = pybamm.ScipySolver(rtol=1e-8, atol=1e-8)
+        t_eval = np.linspace(0, 5, 100)
+        solution = solver.solve(model, t_eval)
+        np.testing.assert_array_almost_equal(
+            solution.y[0], 1 * np.exp(-solution.t), decimal=5
+        )
+
+        # Change initial conditions and solve again
+        model.concatenated_initial_conditions = pybamm.NumpyConcatenation(
+            pybamm.Vector([[2]])
+        )
+        solution = solver.solve(model, t_eval)
+        np.testing.assert_array_almost_equal(
+            solution.y[0], 2 * np.exp(-solution.t), decimal=5
+        )
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")

pip install -e .[dev]
parse_log_pytest
requirements.txt
pytest
apt-get update apt-get install -y gfortran gcc libopenblas-dev liblapack-dev graphviz
3.9
requirements.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@55850464c1b348f61aba25dc38d8c528ae002a6a#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_discretisations/test_discretisation.py::TestDiscretise::test_adding_1D_external_variable tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_concatenation_2D tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_discretise_slicing tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_process_symbol_base tests/unit/test_expression_tree/test_array.py::TestArray::test_list_entries tests/unit/test_expression_tree/test_binary_operators.py::TestIsZero::test_is_matrix_zero tests/unit/test_expression_tree/test_concatenations.py::TestConcatenations::test_base_concatenation tests/unit/test_expression_tree/test_matrix.py::TestMatrix::test_list_entry tests/unit/test_expression_tree/test_operations/test_copy.py::TestCopy::test_symbol_new_copy tests/unit/test_expression_tree/test_operations/test_evaluate.py::TestEvaluate::test_evaluator_python tests/unit/test_expression_tree/test_operations/test_evaluate.py::TestEvaluate::test_find_symbols tests/unit/test_expression_tree/test_operations/test_simplify.py::TestSimplify::test_matrix_simplifications tests/unit/test_expression_tree/test_variable.py::TestVariable::test_variable_bounds tests/unit/test_expression_tree/test_vector.py::TestVector::test_list_entries tests/unit/test_solvers/test_algebraic_solver.py::TestAlgebraicSolver::test_model_solver_least_squares tests/unit/test_solvers/test_algebraic_solver.py::TestAlgebraicSolver::test_model_solver_least_squares_with_bounds tests/unit/test_solvers/test_algebraic_solver.py::TestAlgebraicSolver::test_model_solver_minimize tests/unit/test_solvers/test_algebraic_solver.py::TestAlgebraicSolver::test_model_solver_minimize_with_bounds tests/unit/test_solvers/test_casadi_algebraic_solver.py::TestCasadiAlgebraicSolver::test_model_solver_with_bounds tests/unit/test_solvers/test_scipy_solver.py::TestScipySolver::test_model_solver_manually_update_initial_conditions
tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_mass_matrix_inverse tests/unit/test_solvers/test_base_solver.py::TestBaseSolver::test_fail_consistent_initial_conditions tests/unit/test_solvers/test_casadi_algebraic_solver.py::TestCasadiAlgebraicSolver::test_root_find_fail
tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_add_internal_boundary_conditions tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_adding_0D_external_variable tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_adding_0D_external_variable_fail tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_adding_2D_external_variable_fail tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_broadcast tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_broadcast_2D tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_check_tab_bcs_error tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_concatenate_in_order tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_concatenation tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_concatenation_external_variables tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_concatenation_of_scalars tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_discretise_spatial_operator tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_exceptions tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_no_mesh tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_process_complex_expression tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_process_dict tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_process_empty_model tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_process_input_variable tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_process_model_concatenation tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_process_model_dae tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_process_model_fail tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_process_model_not_inplace tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_process_model_ode tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_process_variables_dict tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_process_with_no_check tests/unit/test_discretisations/test_discretisation.py::TestDiscretise::test_secondary_broadcast_2D tests/unit/test_expression_tree/test_array.py::TestArray::test_linspace tests/unit/test_expression_tree/test_array.py::TestArray::test_meshgrid tests/unit/test_expression_tree/test_array.py::TestArray::test_name tests/unit/test_expression_tree/test_binary_operators.py::TestBinaryOperators::test_addition tests/unit/test_expression_tree/test_binary_operators.py::TestBinaryOperators::test_addition_printing tests/unit/test_expression_tree/test_binary_operators.py::TestBinaryOperators::test_binary_operator tests/unit/test_expression_tree/test_binary_operators.py::TestBinaryOperators::test_binary_operator_domains tests/unit/test_expression_tree/test_binary_operators.py::TestBinaryOperators::test_diff tests/unit/test_expression_tree/test_binary_operators.py::TestBinaryOperators::test_heaviside tests/unit/test_expression_tree/test_binary_operators.py::TestBinaryOperators::test_id tests/unit/test_expression_tree/test_binary_operators.py::TestBinaryOperators::test_inner tests/unit/test_expression_tree/test_binary_operators.py::TestBinaryOperators::test_known_eval tests/unit/test_expression_tree/test_binary_operators.py::TestBinaryOperators::test_minimum_maximum tests/unit/test_expression_tree/test_binary_operators.py::TestBinaryOperators::test_number_overloading tests/unit/test_expression_tree/test_binary_operators.py::TestBinaryOperators::test_power tests/unit/test_expression_tree/test_binary_operators.py::TestBinaryOperators::test_source tests/unit/test_expression_tree/test_binary_operators.py::TestBinaryOperators::test_source_error tests/unit/test_expression_tree/test_binary_operators.py::TestBinaryOperators::test_sparse_divide tests/unit/test_expression_tree/test_binary_operators.py::TestBinaryOperators::test_sparse_multiply tests/unit/test_expression_tree/test_binary_operators.py::TestIsZero::test_is_scalar_zero tests/unit/test_expression_tree/test_concatenations.py::TestConcatenations::test_broadcast_and_concatenate tests/unit/test_expression_tree/test_concatenations.py::TestConcatenations::test_concatenation_auxiliary_domains tests/unit/test_expression_tree/test_concatenations.py::TestConcatenations::test_concatenation_domains tests/unit/test_expression_tree/test_concatenations.py::TestConcatenations::test_concatenation_orphans tests/unit/test_expression_tree/test_concatenations.py::TestConcatenations::test_domain_concatenation_domains tests/unit/test_expression_tree/test_concatenations.py::TestConcatenations::test_domain_error tests/unit/test_expression_tree/test_concatenations.py::TestConcatenations::test_numpy_concatenation_simplify tests/unit/test_expression_tree/test_concatenations.py::TestConcatenations::test_numpy_concatenation_vector_scalar tests/unit/test_expression_tree/test_concatenations.py::TestConcatenations::test_numpy_concatenation_vectors tests/unit/test_expression_tree/test_concatenations.py::TestConcatenations::test_numpy_domain_concatenation tests/unit/test_expression_tree/test_matrix.py::TestMatrix::test_array_wrapper tests/unit/test_expression_tree/test_matrix.py::TestMatrix::test_matrix_evaluate tests/unit/test_expression_tree/test_matrix.py::TestMatrix::test_matrix_operations tests/unit/test_expression_tree/test_operations/test_evaluate.py::TestEvaluate::test_domain_concatenation tests/unit/test_expression_tree/test_operations/test_evaluate.py::TestEvaluate::test_domain_concatenation_2D tests/unit/test_expression_tree/test_operations/test_evaluate.py::TestEvaluate::test_to_python tests/unit/test_expression_tree/test_operations/test_jac_2D.py::TestJacobian::test_functions tests/unit/test_expression_tree/test_operations/test_jac_2D.py::TestJacobian::test_jac_of_domain_concatenation tests/unit/test_expression_tree/test_operations/test_jac_2D.py::TestJacobian::test_linear tests/unit/test_expression_tree/test_operations/test_jac_2D.py::TestJacobian::test_nonlinear tests/unit/test_expression_tree/test_operations/test_simplify.py::TestSimplify::test_domain_concatenation_simplify tests/unit/test_expression_tree/test_operations/test_simplify.py::TestSimplify::test_matrix_divide_simplify tests/unit/test_expression_tree/test_operations/test_simplify.py::TestSimplify::test_simplify_broadcast tests/unit/test_expression_tree/test_operations/test_simplify.py::TestSimplify::test_simplify_concatenation_state_vectors tests/unit/test_expression_tree/test_operations/test_simplify.py::TestSimplify::test_simplify_heaviside tests/unit/test_expression_tree/test_operations/test_simplify.py::TestSimplify::test_simplify_inner tests/unit/test_expression_tree/test_operations/test_simplify.py::TestSimplify::test_symbol_simplify tests/unit/test_expression_tree/test_operations/test_simplify.py::TestSimplify::test_vector_zero_simplify tests/unit/test_expression_tree/test_variable.py::TestVariable::test_variable_diff tests/unit/test_expression_tree/test_variable.py::TestVariable::test_variable_id tests/unit/test_expression_tree/test_variable.py::TestVariable::test_variable_init tests/unit/test_expression_tree/test_variable.py::TestVariableDot::test_variable_diff tests/unit/test_expression_tree/test_variable.py::TestVariableDot::test_variable_id tests/unit/test_expression_tree/test_variable.py::TestVariableDot::test_variable_init tests/unit/test_expression_tree/test_variable.py::TestExternalVariable::test_external_variable_diff tests/unit/test_expression_tree/test_variable.py::TestExternalVariable::test_external_variable_scalar tests/unit/test_expression_tree/test_variable.py::TestExternalVariable::test_external_variable_vector tests/unit/test_expression_tree/test_vector.py::TestVector::test_array_wrapper tests/unit/test_expression_tree/test_vector.py::TestVector::test_column_reshape tests/unit/test_expression_tree/test_vector.py::TestVector::test_vector_evaluate tests/unit/test_expression_tree/test_vector.py::TestVector::test_vector_operations tests/unit/test_expression_tree/test_vector.py::TestVector::test_wrong_size_entries tests/unit/test_solvers/test_algebraic_solver.py::TestAlgebraicSolver::test_algebraic_solver_init tests/unit/test_solvers/test_algebraic_solver.py::TestAlgebraicSolver::test_model_solver tests/unit/test_solvers/test_algebraic_solver.py::TestAlgebraicSolver::test_model_solver_with_time tests/unit/test_solvers/test_algebraic_solver.py::TestAlgebraicSolver::test_root_find_fail tests/unit/test_solvers/test_algebraic_solver.py::TestAlgebraicSolver::test_simple_root_find tests/unit/test_solvers/test_algebraic_solver.py::TestAlgebraicSolver::test_solve_with_input tests/unit/test_solvers/test_algebraic_solver.py::TestAlgebraicSolver::test_with_jacobian tests/unit/test_solvers/test_algebraic_solver.py::TestAlgebraicSolver::test_wrong_solver tests/unit/test_solvers/test_base_solver.py::TestBaseSolver::test_base_solver_init tests/unit/test_solvers/test_base_solver.py::TestBaseSolver::test_block_symbolic_inputs tests/unit/test_solvers/test_base_solver.py::TestBaseSolver::test_convert_to_casadi_format tests/unit/test_solvers/test_base_solver.py::TestBaseSolver::test_discretise_model tests/unit/test_solvers/test_base_solver.py::TestBaseSolver::test_find_consistent_initial_conditions tests/unit/test_solvers/test_base_solver.py::TestBaseSolver::test_nonmonotonic_teval tests/unit/test_solvers/test_base_solver.py::TestBaseSolver::test_ode_solver_fail_with_dae tests/unit/test_solvers/test_base_solver.py::TestBaseSolver::test_root_method_init tests/unit/test_solvers/test_base_solver.py::TestBaseSolver::test_solution_time_length_fail tests/unit/test_solvers/test_base_solver.py::TestBaseSolver::test_step_or_solve_empty_model tests/unit/test_solvers/test_base_solver.py::TestBaseSolver::test_t_eval_none tests/unit/test_solvers/test_base_solver.py::TestBaseSolver::test_timescale_input_fail tests/unit/test_solvers/test_casadi_algebraic_solver.py::TestCasadiAlgebraicSolver::test_algebraic_solver_init tests/unit/test_solvers/test_casadi_algebraic_solver.py::TestCasadiAlgebraicSolver::test_model_solver_with_time tests/unit/test_solvers/test_casadi_algebraic_solver.py::TestCasadiAlgebraicSolver::test_simple_root_find tests/unit/test_solvers/test_casadi_algebraic_solver.py::TestCasadiAlgebraicSolver::test_simple_root_find_correct_initial_guess tests/unit/test_solvers/test_casadi_algebraic_solver.py::TestCasadiAlgebraicSolver::test_solve_with_input tests/unit/test_solvers/test_casadi_algebraic_solver.py::TestCasadiAlgebraicSolverSensitivity::test_least_squares_fit tests/unit/test_solvers/test_casadi_algebraic_solver.py::TestCasadiAlgebraicSolverSensitivity::test_least_squares_fit_input_in_initial_conditions tests/unit/test_solvers/test_casadi_algebraic_solver.py::TestCasadiAlgebraicSolverSensitivity::test_solve_with_symbolic_input tests/unit/test_solvers/test_casadi_algebraic_solver.py::TestCasadiAlgebraicSolverSensitivity::test_solve_with_symbolic_input_1D_scalar_input tests/unit/test_solvers/test_casadi_algebraic_solver.py::TestCasadiAlgebraicSolverSensitivity::test_solve_with_symbolic_input_1D_vector_input tests/unit/test_solvers/test_casadi_algebraic_solver.py::TestCasadiAlgebraicSolverSensitivity::test_solve_with_symbolic_input_in_initial_conditions tests/unit/test_solvers/test_scipy_solver.py::TestScipySolver::test_model_solver_failure tests/unit/test_solvers/test_scipy_solver.py::TestScipySolver::test_model_solver_inputs_in_initial_conditions tests/unit/test_solvers/test_scipy_solver.py::TestScipySolver::test_model_solver_ode_with_jacobian_python tests/unit/test_solvers/test_scipy_solver.py::TestScipySolver::test_model_solver_python tests/unit/test_solvers/test_scipy_solver.py::TestScipySolver::test_model_solver_with_event_python tests/unit/test_solvers/test_scipy_solver.py::TestScipySolver::test_model_solver_with_event_with_casadi tests/unit/test_solvers/test_scipy_solver.py::TestScipySolver::test_model_solver_with_external tests/unit/test_solvers/test_scipy_solver.py::TestScipySolver::test_model_solver_with_inputs tests/unit/test_solvers/test_scipy_solver.py::TestScipySolver::test_model_solver_with_inputs_with_casadi tests/unit/test_solvers/test_scipy_solver.py::TestScipySolver::test_model_step_python
pybamm/expression_tree/array.py pybamm/expression_tree/vector.py pybamm/expression_tree/concatenations.py pybamm/expression_tree/state_vector.py pybamm/spatial_methods/finite_volume.py pybamm/discretisations/discretisation.py pybamm/expression_tree/variable.py pybamm/expression_tree/matrix.py CHANGELOG.md pybamm/solvers/casadi_algebraic_solver.py pybamm/solvers/algebraic_solver.py pybamm/solvers/base_solver.py pybamm/models/standard_variables.py
pybamm/expression_tree/array.py pybamm/expression_tree/vector.py pybamm/expression_tree/concatenations.py pybamm/expression_tree/state_vector.py pybamm/spatial_methods/finite_volume.py pybamm/discretisations/discretisation.py pybamm/expression_tree/variable.py pybamm/expression_tree/matrix.py CHANGELOG.md pybamm/solvers/casadi_algebraic_solver.py pybamm/solvers/algebraic_solver.py pybamm/solvers/base_solver.py pybamm/models/standard_variables.py