diff --git a/process/core/input.py b/process/core/input.py
index 93886b8a16..af1f485909 100644
--- a/process/core/input.py
+++ b/process/core/input.py
@@ -1133,24 +1133,15 @@ def bounds(self) -> tuple[NumberType | None, NumberType | None]:
         "scan",
         int,
         choices=range(IPNSCNS + 1),
+        array=True,
     ),
     "nsweep": InputVariable(
         "scan",
         int,
         choices=range(1, IPNSCNV + 1),
-    ),
-    "isweep_2": InputVariable(
-        "scan",
-        int,
-        choices=range(IPNSCNS + 1),
-    ),
-    "nsweep_2": InputVariable(
-        "scan",
-        int,
-        choices=range(1, IPNSCNV + 1),
+        array=True,
     ),
     "sweep": InputVariable("scan", float, array=True),
-    "sweep_2": InputVariable("scan", float, array=True),
     "impvardiv": InputVariable(
         "reinke",
         int,
diff --git a/process/core/scan.py b/process/core/scan.py
index d62a84a769..cca39075de 100644
--- a/process/core/scan.py
+++ b/process/core/scan.py
@@ -18,11 +18,11 @@
 from process.core.log import logging_model_handler, show_errors
 from process.core.solver import constraints
 from process.core.solver.solver_handler import SolverHandler
-from process.data_structure.scan_variables import IPNSCNS, NOUTVARS, ScanData
 from process.models.availability import AvailabilityModel
 
 if TYPE_CHECKING:
-    from process.core.model import DataStructure, Model
+    from process.core.model import DataStructure
+    from process.main import Models
 
 
 logger = logging.getLogger(__name__)
@@ -224,30 +224,33 @@ def get_val(self, mfile, scan):
     n_tf_coil_turns = (81, Area.T)
 
 
+@dataclass
+class ScanRes:
+    iscan: int
+    ifail: int
+    solver: SolverHandler
+
+
 class Scan:
-    """Perform a parameter scan using the Fortran scan module."""
-
-    def __init__(self, models: Model, solver: str, data: DataStructure):
-        """Immediately run the run_scan() method.
-
-        Parameters
-        ----------
-        models :
-            Physics and engineering model objects
-        solver :
-            Which solver to use, as specified in solver.py
-        data :
-            Data structure object
-        """
+    """Perform a parameter scan
+
+    Parameters
+    ----------
+    models :
+        Physics and engineering model objects
+    solver :
+        Which solver to use, as specified in solver.py
+    data :
+        Data structure object
+    """
+
+    def __init__(self, models: Models, solver: str, data: DataStructure):
         self.models = models
         self.solver = solver
         self.data = data
-        self.solver_handler = SolverHandler(models, solver, data)
-        self.run_scan()
-
-    def run_scan(self):
-        """Call a solver over a range of values of one of the variables.
 
+    def _run(self, iscan, nsweep, sweep, data):
+        """
         This method calls the optimisation routine VMCON a number of times, by
         performing a sweep over a range of values of a particular variable. A
         number of output variable values are written to the MFILE.DAT file at
@@ -258,304 +261,210 @@ def run_scan(self):
         ProcessValueError
             isweep value greater than IPNSCNS
         """
-        if self.data.scan.isweep == 0:
-            # Solve single problem, rather than an array of problems (scan)
-            # doopt() can also run just an evaluation
-            start_time = time.time()
-            ifail = self.doopt()
-            write_output_files(
-                models=self.models,
-                data=self.data,
-                ifail=ifail,
-                runtime=time.time() - start_time,
-            )
-            show_errors(constants.NOUT)
-            return
-
-        if self.data.scan.isweep > IPNSCNS:
-            raise ProcessValueError(
-                "Illegal value of isweep",
-                isweep=self.data.scan.isweep,
-                IPNSCNS=IPNSCNS,
-            )
+        sh = SolverHandler(self.models, self.solver, data)
+        # TODO queue the output to avoid race condition (?)
+        if data.scan.nsweep is not None:
+            self.write_point_header(iscan)
+        start_time = time.time()
+        ifail = sh.run()
+        end_time = time.time() - start_time
+        write_output_files(models=self.models, data=data, ifail=ifail, runtime=end_time)
+        nums = data.numerics
+        nums.sqsumsq = sum(r**2 for r in nums.rcm[: nums.neqns]) ** 0.5
 
         if self.data.scan.scan_dim == 2:
             self.scan_2d()
         else:
             self.scan_1d()
 
-    def doopt(self):
-        """Run the optimiser or solver."""
-        ifail = self.solver_handler.run()
-        constraints.constraints_output(self.data, self.solver)
-
-        return ifail
-
-    def scan_1d(self):
-        """Run a 1-D scan."""
-        # initialise dict which will contain ifail values for each scan point
-        scan_1d_ifail_dict = {}
-
-        for iscan in range(1, self.data.scan.isweep + 1):
-            self.scan_1d_write_point_header(iscan)
-            start_time = time.time()
-            ifail = self.doopt()
-            scan_1d_ifail_dict[iscan] = ifail
-            write_output_files(
-                models=self.models,
-                data=self.data,
-                ifail=ifail,
-                runtime=time.time() - start_time,
-            )
+    def write_outputs(self):
+        write_output_files(
+            models=self.models,
+            data=self.data,
+            ifail=self._ifail,
+            runtime=self._finish_time - self._start_time,
+        )
+        show_errors(constants.NOUT)
 
-            show_errors(constants.NOUT)
-            logging_model_handler.clear_logs()
+        logging_model_handler.clear_logs()
+        optimisation_output(data)
+        constraints.constraints_output(data, self.solver)
 
-        # outvar now contains results
-        self.scan_1d_write_plot(self.data.scan)
-        print("Scan Convergence Summary \n")
-        sweep_values = self.data.scan.sweep[: self.data.scan.isweep]
-        nsweep_var = self.scan_select(
-            self.data.scan.nsweep, self.data.scan.sweep, self.data.scan.isweep
-        )
-        converged_count = 0
-        # offsets for aligning the converged/unconverged column
-        max_sweep_value_length = len(str(np.max(sweep_values)).replace(".", ""))
-        offsets = [
-            max_sweep_value_length - len(str(sweep_val).replace(".", ""))
-            for sweep_val in sweep_values
+        return ScanRes(iscan, ifail, sh)
+
+    def _set_v_x_label(self, iscan: list[int]):
+        sv = [
+            self.scan_select(self.data.scan.nsweep, self.data.scan.sweep, isc)
+            for isc in iscan
         ]
-        for iscan in range(1, self.data.scan.isweep + 1):
-            if scan_1d_ifail_dict[iscan] == 1:
-                converged_count += 1
-                print(
-                    f"Scan {iscan:02d}: {nsweep_var.fname} = {sweep_values[iscan - 1]} "
-                    + " " * offsets[iscan - 1]
-                    + "\u001b[32mCONVERGED \u001b[0m"
-                )
-            else:
-                print(
-                    f"Scan {iscan:02d}: {nsweep_var.fname} = {sweep_values[iscan - 1]} "
-                    + " " * offsets[iscan - 1]
-                    + "\u001b[31mUNCONVERGED \u001b[0m"
-                )
-        converged_percentage = converged_count / self.data.scan.isweep * 100
-        print(f"\nConvergence Percentage: {converged_percentage:.2f}%")
+        self.data.globals.vlabel = [s.fname for s in sv]
+        self.data.globals.xlabel = [s.data.description for s in sv]
 
-    def scan_2d(self):
-        """Run a 2-D scan."""
-        # Initialise intent(out) arrays
-        self.scan_2d_init(self.data.scan)
-        iscan = 1
-
-        # initialise array which will contain ifail values for each scan point
-        scan_2d_ifail_list = np.zeros(
-            (NOUTVARS, IPNSCNS),
-            dtype=np.float64,
-            order="F",
-        )
-        for iscan_1 in range(1, self.data.scan.isweep + 1):
-            for iscan_2 in range(1, self.data.scan.isweep_2 + 1):
-                self.scan_2d_write_point_header(iscan, iscan_1, iscan_2)
-                start_time = time.time()
-                ifail = self.doopt()
-                write_output_files(
-                    models=self.models,
-                    data=self.data,
-                    ifail=ifail,
-                    runtime=time.time() - start_time,
-                )
+    def write_point_header(self, iscan):
+        self._set_v_x_label(iscan)
 
-                show_errors(constants.NOUT)
-                logging_model_handler.clear_logs()
-                scan_2d_ifail_list[iscan_1][iscan_2] = ifail
-                iscan += 1
+        process_output.oblnkl(constants.NOUT)
+        process_output.oblnkl(constants.MFILE)
 
-        print("Scan Convergence Summary\n")
-        sweep_1_values = self.data.scan.sweep[: self.data.scan.isweep]
-        sweep_2_values = self.data.scan.sweep_2[: self.data.scan.isweep_2]
-        nsweep_var = self.scan_select(
-            self.data.scan.nsweep, self.data.scan.sweep, self.data.scan.isweep
-        )
-        nsweep_2_var = self.scan_select(
-            self.data.scan.nsweep_2, self.data.scan.sweep_2, self.data.scan.isweep_2
-        )
-        converged_count = 0
-        scan_point = 1
-        # offsets for aligning the converged/unconverged column
-        max_sweep1_value_length = len(str(np.max(sweep_1_values)).replace(".", ""))
-        max_sweep2_value_length = len(str(np.max(sweep_2_values)).replace(".", ""))
-        offsets = np.zeros(
-            (self.data.scan.isweep, self.data.scan.isweep_2), dtype=int, order="F"
+        process_output.write(
+            constants.NOUT,
+            f"Scan point {iscan} of {np.prod(self.data.scan.isweep)} : \n".join(
+                f"{v} = {self.data.scan.sweep[iscan[no] - 1]}"
+                for no, v in enumerate(self.data.globals.vlabel)
+            ),
         )
-        for count1, sweep1 in enumerate(sweep_1_values):
-            for count2, sweep2 in enumerate(sweep_2_values):
-                offsets[count1][count2] = (
-                    max_sweep1_value_length
-                    - len(str(sweep1).replace(".", ""))
-                    + max_sweep2_value_length
-                    - len(str(sweep2).replace(".", ""))
-                )
+        process_output.ovarin(constants.MFILE, "Scan point number", "(iscan)", iscan)
 
-        for iscan_1 in range(1, self.data.scan.isweep + 1):
-            for iscan_2 in range(1, self.data.scan.isweep_2 + 1):
-                if scan_2d_ifail_list[iscan_1][iscan_2] == 1:
-                    converged_count += 1
-                    print(
-                        (
-                            f"Scan {scan_point:02d}: ({nsweep_var.fname} = "
-                            f"{sweep_1_values[iscan_1 - 1]}, {nsweep_2_var.fname} "
-                            f"= {sweep_2_values[iscan_2 - 1]}) "
-                        )
-                        + " " * offsets[iscan_1 - 1][iscan_2 - 1]
-                        + "\u001b[32mCONVERGED \u001b[0m"
-                    )
-                    scan_point += 1
-                else:
-                    print(
-                        (
-                            f"Scan {scan_point:02d}: ({nsweep_var.fname} = "
-                            f"{sweep_1_values[iscan_1 - 1]}, {nsweep_2_var.fname} = "
-                            f"{sweep_2_values[iscan_2 - 1]}) "
-                        )
-                        + " " * offsets[iscan_1 - 1][iscan_2 - 1]
-                        + "\u001b[31mUNCONVERGED \u001b[0m"
-                    )
-                    scan_point += 1
-        converged_percentage = (
-            converged_count / (self.data.scan.isweep * self.data.scan.isweep_2) * 100
+        print(
+            f"Starting scan point {iscan}:  {self.data.globals.xlabel}, \n".join(
+                f"{v} = {self.data.scan.sweep[iscan[no] - 1]}"
+                for no, v in enumerate(self.data.globals.vlabel)
+            )
         )
-        print(f"\nConvergence Percentage: {converged_percentage:.2f}%")
 
-    @staticmethod
-    def scan_2d_init(scan_data: ScanData):
-        """Scan 2d initialisation"""
-        process_output.ovarre(
-            constants.MFILE,
-            "Number of first variable scan points",
-            "(isweep)",
-            scan_data.isweep,
-        )
-        process_output.ovarre(
-            constants.MFILE,
-            "Number of second variable scan points",
-            "(isweep_2)",
-            scan_data.isweep_2,
-        )
-        process_output.ovarre(
-            constants.MFILE,
-            "Scanning first variable number",
-            "(nsweep)",
-            scan_data.nsweep,
-        )
-        process_output.ovarre(
-            constants.MFILE,
-            "Scanning second variable number",
-            "(nsweep_2)",
-            scan_data.nsweep_2,
-        )
-        process_output.ovarre(
-            constants.MFILE,
-            "Scanning second variable number",
-            "(nsweep_2)",
-            scan_data.nsweep_2,
-        )
-        process_output.ovarre(
-            constants.MFILE,
-            "Scanning second variable number",
-            "(nsweep_2)",
-            scan_data.nsweep_2,
-        )
+    def scan_select(self, nsweep, sweep, iscan):
+        sv = ScanVariables(nsweep)
+        sv.set(self.data, sweep[iscan - 1])
+        return sv
 
-    def scan_1d_write_point_header(self, iscan: int):
-        """Scan 1d header"""
-        self.data.globals.iscan_global = iscan
-        sv = self.scan_select(self.data.scan.nsweep, self.data.scan.sweep, iscan)
+    def run(self):
+        """Call a solver over a range of values of one of the variables.
+
+        This method calls the optimisation routine VMCON a number of times, by
+        performing a sweep over a range of values of a particular variable. A
+        number of output variable values are written to the MFILE.DAT file at
+        each scan point, for plotting or other post-processing purposes.
+        """
+        # vectorise running of self._run
+        if self.data.scan.nsweep is not None:
+            for d, n, v in (
+                ("Number of scan points", "(isweep)", self.data.scan.isweep),
+                ("Scanning variable number", "(nsweep)", self.data.scan.nsweep),
+            ):
+                process_output.ovarin(constants.MFILE, d, n, v)
+
+        # TODO copy of self.data for each vectorised run (?)
+        scan_res = np.vectorise(self._run)(
+            self.data.scan.isweep, self.data.scan.nsweep, self.data.scan.sweep, self.data
+        )
 
-        self.data.globals.vlabel = sv.fname
-        self.data.globals.xlabel = sv.description
+        if self.data.scan.nsweep is not None:
+            self.summary(scan_res)
 
-        process_output.oblnkl(constants.NOUT)
-        process_output.ostars(constants.NOUT, 110)
+    def summary(self, scan_res):
+        print("Scan Convergence Summary\n")
+        sweep_values = self.data.scan.sweep
+        nsweep_var = [ScanVariables(nsw) for nsw in self.data.scan.nsweep]
 
-        process_output.write(
-            constants.NOUT,
-            f"***** Scan point {iscan} of {self.data.scan.isweep} : "
-            f"{self.data.globals.xlabel}"
-            f", {self.data.globals.vlabel} = {self.data.scan.sweep[iscan - 1]} "
-            "*****",
-        )
-        process_output.ostars(constants.NOUT, 110)
-        process_output.oblnkl(constants.MFILE)
-        process_output.ovarre(constants.MFILE, "Scan point number", "(iscan)", iscan)
+        conv_list = []
+        converged_count = 0
+        conv_str = "\u001b[3{}CONVERGED \u001b[0m"
+        for no, sr in enumerate(scan_res):
+            if sr.ifail == 1:
+                converged_count += 1
+                conv = conv_str.format("2m")
+            else:
+                conv = conv_str.format("1mUN")
+            conv_list.append([
+                "{sr.iscan:02d}",
+                nsweep_var[no].fname,
+                sweep_values[sr.iscan],
+                conv,
+            ])
 
         print(
-            f"Starting scan point {iscan} of {self.data.scan.isweep} : "
-            f"{self.data.globals.xlabel} , {self.data.globals.vlabel}"
-            f" = {self.data.scan.sweep[iscan - 1]}"
+            tabulate(conv_list, headers=["Iscan", "Sweep Var", "Sweep Val", "Converged"])
         )
 
-    def scan_2d_write_point_header(self, iscan, iscan_1, iscan_2):
-        """Scan 2d header"""
-        iscan_r = self.data.scan.isweep_2 - iscan_2 + 1 if iscan_1 % 2 == 0 else iscan_2
+        converged_percentage = converged_count / np.prod(self.data.scan.isweep) * 100
+        print(f"\nConvergence Percentage: {converged_percentage:.2f}%")
 
-        # Makes iscan available globally (read-only)
-        self.data.globals.iscan_global = iscan
-        sv_1 = self.scan_select(self.data.scan.nsweep, self.data.scan.sweep, iscan_1)
 
-        self.data.globals.vlabel = sv_1.fname
-        self.data.globals.xlabel = sv_1.data.description
+def optimisation_output(data: DataStructure):
+    nums = data.numerics
 
-        sv_2 = self.scan_select(self.data.scan.nsweep_2, self.data.scan.sweep_2, iscan_r)
+    written_warning = False
 
-        self.data.globals.vlabel_2 = sv_2.fname
-        self.data.globals.xlabel_2 = sv_2.data.description
+    # Output optimisation parameters
+    solution_vector_table = []
+    for i in range(nums.nvar):
+        nums.xcs[i] = nums.xcm[i] * nums.scafc[i]
 
-        process_output.oblnkl(constants.NOUT)
-        process_output.ostars(constants.NOUT, 110)
+        name = nums.lablxc[nums.ixc[i] - 1]
+        solution_vector_table.append([name, nums.xcs[i], nums.xcm[i]])
 
-        process_output.write(
-            constants.NOUT,
-            f"***** 2D Scan point {iscan} of "
-            f"{self.data.scan.isweep * self.data.scan.isweep_2} : "
-            f"{self.data.globals.vlabel} = {self.data.scan.sweep[iscan_1 - 1]} and"
-            f" {self.data.globals.vlabel_2} = {self.data.scan.sweep_2[iscan_r - 1]} "
-            "*****",
-        )
-        process_output.ostars(constants.NOUT, 110)
-        process_output.oblnkl(constants.MFILE)
-        process_output.ovarre(constants.MFILE, "Scan point number", "(iscan)", iscan)
+        xminn = 1.01 * nums.itv_scaled_lower_bounds[i]
+        xmaxx = 0.99 * nums.itv_scaled_upper_bounds[i]
 
-        print(
-            f"Starting scan point {iscan}:  {self.data.globals.xlabel}, "
-            f"{self.data.globals.vlabel} = {self.data.scan.sweep[iscan_1 - 1]}"
-            f" and {self.data.globals.xlabel_2}, "
-            f"{self.data.globals.vlabel_2} = {self.data.scan.sweep_2[iscan_r - 1]} "
-        )
+        # Write to output file if close to optimisation parameter bounds
+        if nums.xcm[i] < xminn or nums.xcm[i] > xmaxx:
+            if not written_warning:
+                written_warning = True
+                process_output.ocmmnt(
+                    constants.NOUT,
+                    (
+                        "Certain operating limits have been reached,"
+                        "\n as shown by the following optimisation parameters that are"
+                        "\n at or near to the edge of their prescribed range :\n"
+                    ),
+                )
 
-        return iscan_r
-
-    @staticmethod
-    def scan_1d_write_plot(scan_data: ScanData):
-        """Scan 1d plotter"""
-        if scan_data.first_call_1d:
-            process_output.ovarre(
-                constants.MFILE,
-                "Number of scan points",
-                "(isweep)",
-                scan_data.isweep,
-            )
-            process_output.ovarre(
-                constants.MFILE,
-                "Scanning variable number",
-                "(nsweep)",
-                scan_data.nsweep,
+            xcval = nums.xcm[i] * nums.scafc[i]
+
+            if nums.xcm[i] < xminn:
+                location, bound = "below", "lower"
+                bounds = nums.itv_scaled_lower_bounds
+            else:
+                location, bound = "above", "upper"
+                bounds = nums.itv_scaled_upper_bounds
+            process_output.write(
+                constants.NOUT,
+                f"   {name:<30}= {xcval} is at or {location} its {bound} bound:"
+                f" {bounds[i] * nums.scafc[i]}",
             )
 
-            scan_data.first_call_1d = False
+        xnorm = (
+            1.0
+            if nums.boundu[i] == nums.boundl[i]
+            else min(
+                max(
+                    (nums.xcm[i] - nums.itv_scaled_lower_bounds[i])
+                    / (
+                        nums.itv_scaled_upper_bounds[i] - nums.itv_scaled_lower_bounds[i]
+                    ),
+                    0.0,
+                ),
+                1.0,
+            )
+        )
 
-    def scan_select(self, nsweep, sweep, iscan) -> ScanVariables:
-        """Select a scan"""
-        sv = ScanVariables(nsweep)
-        sv.set(self.data, sweep[iscan - 1])
-        return sv
+        # Write optimisation parameters to mfile
+        for d, var, v in (
+            (nums.lablxc[nums.ixc[i] - 1], f"(itvar{i + 1:03d})", nums.xcs[i]),
+            (f"{name} (final value/initial value)", f"(xcm{i + 1:03d})", nums.xcm[i]),
+            (f"{name} (range normalised)", f"(nitvar{i + 1:03d})", xnorm),
+            (
+                f"{name} (upper bound)",
+                f"(boundu{i + 1:03d})",
+                nums.itv_scaled_upper_bounds[i] * nums.scafc[i],
+            ),
+            (
+                f"{name} (lower bound)",
+                f"(boundl{i + 1:03d})",
+                nums.itv_scaled_lower_bounds[i] * nums.scafc[i],
+            ),
+        ):
+            process_output.ovarre(constants.MFILE, d, var, v)
+
+    # Write optimisation parameter headings to output file
+    process_output.osubhd(
+        constants.NOUT, "The solution vector is comprised as follows :"
+    )
+    process_output.write(
+        constants.NOUT,
+        tabulate(
+            solution_vector_table,
+            headers=["", "Final value", "Final / initial"],
+            numalign="left",
+        ),
+    )
diff --git a/process/core/solver/iteration_variables.py b/process/core/solver/iteration_variables.py
index c40f932ec1..85b0c07d4c 100644
--- a/process/core/solver/iteration_variables.py
+++ b/process/core/solver/iteration_variables.py
@@ -326,10 +326,7 @@ def load_iteration_variables(data):
 
         # warn of the iteration variable is also a scan variable because this will cause
         # the optimiser and scan to overwrite the same variable and conflict
-        if iteration_variable.name in {
-            data.globals.vlabel,
-            data.globals.vlabel_2,
-        }:
+        if iteration_variable.name in data.globals.vlabel:
             logger.critical(
                 (
                     "The sweep variable is also an iteration variable and will be "
diff --git a/process/core/solver/solver.py b/process/core/solver/solver.py
index 269177e338..513bebeafa 100644
--- a/process/core/solver/solver.py
+++ b/process/core/solver/solver.py
@@ -325,11 +325,6 @@ def verror(self):
                     "A feasible solution may be difficult to achieve.",
                     "Try changing or adding variables to IXC.",
                 ),
-                4: (
-                    "An uphill search direction was found.",
-                    "Try changing the equations in ICC, or",
-                    "adding new variables to IXC.",
-                ),
                 SolverOutputCondition.NO_SOLUTION: (
                     "The quadratic programming technique was unable to",
                     "find a feasible point.\n",
@@ -337,13 +332,6 @@ def verror(self):
                     "their initial values (especially if only 1 optimisation",
                     "iteration was performed).",
                 ),
-                6: (
-                    "The quadratic programming technique was restricted",
-                    "by an artificial bound, or failed due to a singular",
-                    "matrix.",
-                    "Try changing the equations in ICC, or",
-                    "adding new variables to IXC.",
-                ),
             }.get(self.info, "Unknown Error code")
         )
 
diff --git a/process/core/solver/solver_handler.py b/process/core/solver/solver_handler.py
index 4901d5983a..6c4cb30c87 100644
--- a/process/core/solver/solver_handler.py
+++ b/process/core/solver/solver_handler.py
@@ -1,6 +1,7 @@
 """Module containing solver handler routines"""
 
 import logging
+from contextlib import contextmanager
 
 from tabulate import tabulate
 
@@ -50,14 +51,9 @@ def run(self):
         # Initialise iteration variables and bounds in Python: relies on Fortran
         # iteration variables being defined above
         # Trim maximum size arrays down to actually used size
-        n = self.data.numerics.nvar
-        x = self.data.numerics.xcm[:n]
-        bndl = self.data.numerics.itv_scaled_lower_bounds[:n]
-        bndu = self.data.numerics.itv_scaled_upper_bounds[:n]
-
-        # Define total number of constraints and equality constraints
-        m = self.data.numerics.neqns + self.data.numerics.nineqns
-        meq = self.data.numerics.neqns
+        x = self.data.numerics.xcm[: self.data.numerics.nvar]
+        bndl = self.data.numerics.itv_scaled_lower_bounds[: self.data.numerics.nvar]
+        bndu = self.data.numerics.itv_scaled_upper_bounds[: self.data.numerics.nvar]
 
         # Evaluators() calculates the objective and constraint functions and
         # their gradients for a given vector x
@@ -68,31 +64,18 @@ def run(self):
         self.solver.set_evaluators(evaluators)
         self.solver.set_bounds(bndl, bndu)
         self.solver.set_opt_params(x)
-        self.solver.set_constraints(m, meq)
+        # Define total number of constraints and equality constraints
+        self.solver.set_constraints(
+            m=self.data.numerics.neqns + self.data.numerics.nineqns,
+            meq=self.data.numerics.neqns,
+        )
         ifail = self.solver.solve()
 
         # If VMCON optimisation has failed then try altering value of epsfcn
         if self.solver_name == "vmcon":
-            if ifail != SolverOutputCondition.CONVERGED:
-                print("Trying again with new epsfcn")
-                # epsfcn is only used in evaluators.Evaluators()
-                # TODO epsfcn could be set in Evaluators instance now, don't need to
-                # set/unset in self.data.numerics module
-                self.data.numerics.epsfcn *= 10  # try new larger value
-                print("new epsfcn = ", self.data.numerics.epsfcn)
-
-                ifail = self.solver.solve()
-                # First solution attempt failed
-                # (ifail != SolverOutputCondition.CONVERGED): supply ifail value
-                # to next attempt
-                self.data.numerics.epsfcn /= 10  # reset value
-
-            if ifail != SolverOutputCondition.CONVERGED:
-                print("Trying again with new epsfcn")
-                self.data.numerics.epsfcn /= 10  # try new smaller value
-                print("new epsfcn = ", self.data.numerics.epsfcn)
-                ifail = self.solver.solve()
-                self.data.numerics.epsfcn *= 10  # reset value
+            if ifail != 1:
+                with epsfcn_context(self.data.numerics):
+                    self.solver.solve()
 
             # If VMCON has exited with error code 5
             # (ifail = SolverOutputCondition.NO_SOLUTION) try another run using a
@@ -346,3 +329,17 @@ def _optimisation_parameters_output(self):
                 numalign="left",
             ),
         )
+
+
+@contextmanager
+def epsfcn_context(numerics):
+    print("Trying again with new epsfcn")
+    # epsfcn is only used in evaluators.Evaluators()
+    # TODO epsfcn could be set in Evaluators instance now, don't need to
+    # set/unset in numerics module
+    numerics.epsfcn *= 10  # try new larger value
+    print("new epsfcn = ", numerics.epsfcn)
+    try:
+        yield
+    finally:
+        numerics.epsfcn /= 10  # reset value
diff --git a/process/data_structure/global_variables.py b/process/data_structure/global_variables.py
index f00c9441ab..a84275727e 100644
--- a/process/data_structure/global_variables.py
+++ b/process/data_structure/global_variables.py
@@ -1,6 +1,6 @@
 """Module containing variables necessary for running PROCESS"""
 
-from dataclasses import dataclass
+from dataclasses import dataclass, field
 
 
 @dataclass(slots=True)
@@ -22,19 +22,11 @@ class GlobalData:
     output_prefix: str = ""
     """Output file path prefix"""
 
-    xlabel: str = ""
-    """Scan parameter description label"""
+    xlabel: list[str] = field(default_factory=lambda: [""])
+    """Scan parameters description label"""
 
-    vlabel: str = ""
-    """Scan value name label"""
-
-    xlabel_2: str = ""
-    """Scan parameter description label (2nd dimension)"""
-
-    vlabel_2: str = ""
-    """Scan value name label (2nd dimension)"""
-
-    iscan_global: int = 0
+    vlabel: list[str] = field(default_factory=lambda: [""])
+    """Scan values name label"""
 
     convergence_parameter: float = 0.0
     """VMCON convergence parameter 'sum'"""
diff --git a/process/data_structure/scan_variables.py b/process/data_structure/scan_variables.py
index 0a268e81e1..18d06ccf4e 100644
--- a/process/data_structure/scan_variables.py
+++ b/process/data_structure/scan_variables.py
@@ -4,17 +4,12 @@
 
 import numpy as np
 
+from process.core.exceptions import ProcessValueError
+
 IPNSCNS = 1000
 """Maximum number of scan points"""
 
 
-IPNSCNV = 81
-"""Number of available scan variables"""
-
-
-NOUTVARS = 84
-
-
 @dataclass(slots=True)
 class ScanData:
     """Dataclass holding scan variables"""
@@ -22,35 +17,45 @@ class ScanData:
     scan_dim: int = 1
     """1-D or 2-D scan switch (1=1D, 2=2D)"""
 
-    isweep: int = 0
+    isweep: list[int] | int = 1
     """Number of scan points to calculate"""
 
-    isweep_2: int = 0
-    """Number of 2D scan points to calculate"""
-
-    nsweep: int = 1
+    nsweep: list[int] | int | None = None
     """Switch denoting quantity to scan
+    see `process.core.scan.ScanVariables` for available options
     """
 
-    nsweep_2: int = 3
-    """nsweep_2 /3/ : switch denoting quantity to scan for 2D scan:"""
-
-    sweep: list[float] = field(
-        default_factory=lambda: np.zeros(IPNSCNS, dtype=np.float64)
-    )
-    """sweep(IPNSCNS) /../: actual values to use in scan"""
-
-    sweep_2: list[float] = field(
-        default_factory=lambda: np.zeros(IPNSCNS, dtype=np.float64)
-    )
-    """sweep_2(IPNSCNS) /../: actual values to use in 2D scan"""
-
-    # Vars in subroutines scan_1d and scan_2d requiring re-initialising before
-    # each new run
-
-    first_call_1d: bool = True
-
-    first_call_2d: bool = True
+    sweep: np.ndarray = field(default_factory=lambda: np.zeros(1, dtype=np.float64))
+    """Actual values to use in scan"""
+
+    def __post_init__(self):
+        if isinstance(self.isweep, int):
+            # avoid old 0 default
+            self.isweep = [self.isweep or 1]
+
+        if len(self.isweep) > 2 or len(self.sweep.shape) > 2:
+            raise NotImplementedError("N-D Scans not currently supported")
+
+        if max(self.isweep) > IPNSCNS:
+            raise ProcessValueError(
+                "Illegal value of isweep",
+                isweep=self.isweep,
+                IPNSCNS=IPNSCNS,
+            )
+        if self.nsweep != len(self.isweep):
+            raise ValueError(
+                "Number of sweep variables not equal to scan point dimensions"
+            )
+        if self.sweep.shape != self.isweep:
+            if self.isweep != 1:
+                self.isweep = list(self.sweep.shape)
+            else:
+                print("Unset sweep values set to zero")
+                # TODO append to size instead of resetting
+                self.sweep = np.zeros(self.isweep, dtype=np.float64)
+
+        self.nsweep = np.asarray(self.nsweep, dtype=int)
+        self.isweep = np.asarray(self.isweep, dtype=int)
 
 
 CREATE_DICTS_FROM_DATACLASS = ScanData
diff --git a/process/main.py b/process/main.py
index 2501505685..2f326b924c 100644
--- a/process/main.py
+++ b/process/main.py
@@ -459,6 +459,7 @@ def run_scan(self):
                 "select either 1 (optimise) or -2 (no optimisation)."
             )
         self.scan = Scan(self.models, self.solver, self.data)
+        self.scan.run()
 
     @staticmethod
     def show_errors():