diff --git a/CHANGELOG.md b/CHANGELOG.md
index 4ffe632a82..a0c480ad92 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 2afdb91020..1e6de70f92 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -160,8 +160,6 @@ def process_model(self, model, inplace=True, check_model=True):
         # 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 @@ def process_model(self, model, inplace=True, check_model=True):
             # 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 @@ def set_variable_slices(self, variables):
         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 @@ def set_variable_slices(self, variables):
                     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 @@ def set_variable_slices(self, variables):
         # 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 @@ def _process_symbol(self, symbol):
                 # 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 @@ def _process_symbol(self, symbol):
             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 @@ def _process_symbol(self, symbol):
                 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 41bed0708d..95622571b4 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 @@ def __init__(
         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 611b45773b..456029314a 100644
--- a/pybamm/expression_tree/concatenations.py
+++ b/pybamm/expression_tree/concatenations.py
@@ -127,7 +127,7 @@ def __init__(self, *children):
         # 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 cd613266de..8aed40bc11 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 @@ def __init__(
         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 47a6c28e11..7bca90e6c1 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 3be2f0b614..44055c651d 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 7714ab1769..a03f177f4f 100644
--- a/pybamm/expression_tree/vector.py
+++ b/pybamm/expression_tree/vector.py
@@ -21,6 +21,8 @@ def __init__(
         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 bf0b995583..87e0c7b03f 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 @@
         "secondary": "negative electrode",
         "tertiary": "current collector",
     },
+    bounds=(0, 1),
 )
 c_s_p = pybamm.Variable(
     "Positive particle concentration",
@@ -113,32 +120,41 @@
         "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 @@
     "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 @@
     domain=["negative electrode"],
     auxiliary_domains={"secondary": "current collector"},
 )
-
diff --git a/pybamm/solvers/algebraic_solver.py b/pybamm/solvers/algebraic_solver.py
index da93e50d9e..dd4e0d0b20 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 @@ def jac_fn(y_alg):
                 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 @@ def jac_fn(y_alg):
                     )
                 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 159118bf72..3217097817 100644
--- a/pybamm/solvers/base_solver.py
+++ b/pybamm/solvers/base_solver.py
@@ -51,7 +51,7 @@ def __init__(
                 "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 @@ def copy(self):
         "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 @@ def solve(self, model, t_eval=None, external_variables=None, inputs=None):
         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 229cc7f6dd..e286936f24 100644
--- a/pybamm/solvers/casadi_algebraic_solver.py
+++ b/pybamm/solvers/casadi_algebraic_solver.py
@@ -68,6 +68,7 @@ def _integrate(self, model, t_eval, inputs=None):
         # 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 @@ def _integrate(self, model, t_eval, inputs=None):
         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 @@ def _integrate(self, model, t_eval, inputs=None):
                         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 c2e9683cbd..8b00a6bfc5 100644
--- a/pybamm/spatial_methods/finite_volume.py
+++ b/pybamm/spatial_methods/finite_volume.py
@@ -641,7 +641,7 @@ def add_ghost_nodes(self, symbol, discretised_symbol, bcs):
         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"]
diff --git a/tests/unit/test_discretisations/test_discretisation.py b/tests/unit/test_discretisations/test_discretisation.py
index f827a5a135..7ceff8746c 100644
--- a/tests/unit/test_discretisations/test_discretisation.py
+++ b/tests/unit/test_discretisations/test_discretisation.py
@@ -131,6 +131,7 @@ def test_adding_1D_external_variable(self):
         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 @@ def test_discretise_slicing(self):
         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 @@ def test_discretise_slicing(self):
             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 @@ def test_discretise_slicing(self):
                 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 @@ def test_process_symbol_base(self):
         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 @@ def test_process_model_dae(self):
         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 @@ def test_concatenation_2D(self):
             "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 @@ def test_concatenation_2D(self):
         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 550abd84c9..e2adab03b7 100644
--- a/tests/unit/test_expression_tree/test_array.py
+++ b/tests/unit/test_expression_tree/test_array.py
@@ -12,6 +12,12 @@ def test_name(self):
         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 02a4b5ceab..71d1689c33 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 @@ def test_is_matrix_zero(self):
 
         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 3e2e42b2b2..e87f2aaa29 100644
--- a/tests/unit/test_expression_tree/test_concatenations.py
+++ b/tests/unit/test_expression_tree/test_concatenations.py
@@ -18,9 +18,9 @@ def test_base_concatenation(self):
         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 0d2300ed53..f4eb8988d9 100644
--- a/tests/unit/test_expression_tree/test_matrix.py
+++ b/tests/unit/test_expression_tree/test_matrix.py
@@ -19,6 +19,12 @@ def test_array_wrapper(self):
         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 ef6d3a0638..564eee5e61 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 @@ def test_symbol_new_copy(self):
         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 8ff7eaaf43..3f08f11374 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 @@ def test_find_symbols(self):
         # 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 @@ def test_evaluator_python(self):
             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 @@ def test_evaluator_python(self):
             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 @@ def test_evaluator_python(self):
             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 @@ def test_evaluator_python(self):
             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 @@ def test_evaluator_python(self):
             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 3bd80db91b..5fea85ee74 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 @@ def test_jac_of_domain_concatenation(self):
         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 b74259d89a..3aac02b6f9 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 @@ def test_matrix_simplifications(self):
         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 @@ def test_matrix_simplifications(self):
         )
 
         # 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 @@ def test_matrix_simplifications(self):
             )
 
         # 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 8e4a8be11f..b023f1e439 100644
--- a/tests/unit/test_expression_tree/test_variable.py
+++ b/tests/unit/test_expression_tree/test_variable.py
@@ -33,6 +33,18 @@ def test_variable_id(self):
         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 924876dd83..52c414b78c 100644
--- a/tests/unit/test_expression_tree/test_vector.py
+++ b/tests/unit/test_expression_tree/test_vector.py
@@ -21,6 +21,12 @@ def test_column_reshape(self):
         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 45c858327c..4e7d660b55 100644
--- a/tests/unit/test_solvers/test_algebraic_solver.py
+++ b/tests/unit/test_solvers/test_algebraic_solver.py
@@ -135,6 +135,7 @@ def test_model_solver(self):
 
         # 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 @@ def test_model_solver(self):
             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 3633cf6bf1..35e6cf25d7 100644
--- a/tests/unit/test_solvers/test_base_solver.py
+++ b/tests/unit/test_solvers/test_base_solver.py
@@ -117,6 +117,7 @@ def __init__(self):
                     "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 @@ def __init__(self):
                     "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 @@ def __init__(self):
                     "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 208860d047..1d8422318f 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 Model:
             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 @@ def test_model_solver_with_time(self):
             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 c984d7f408..6ab8fe2dcb 100644
--- a/tests/unit/test_solvers/test_scikits_solvers.py
+++ b/tests/unit/test_solvers/test_scikits_solvers.py
@@ -42,7 +42,7 @@ def test_model_dae_integrate_failure_bad_ics(self):
 
         # 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 @@ def test_dae_integrate_with_non_unity_mass(self):
 
         # 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 1d6cfd43a6..1d0aae34b5 100644
--- a/tests/unit/test_solvers/test_scipy_solver.py
+++ b/tests/unit/test_solvers/test_scipy_solver.py
@@ -314,6 +314,30 @@ def test_model_solver_inputs_in_initial_conditions(self):
             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")