diff --git a/docs/source/tutorial.rst b/docs/source/tutorial.rst
index 7fd7db272ae61b147c5ea489ed5f380ca6e69ce7..134e2c3d5f5a861d70933a1aec4c423f7d93957e 100644
--- a/docs/source/tutorial.rst
+++ b/docs/source/tutorial.rst
@@ -83,9 +83,8 @@ In python this could be implemented like this :class:`~solver.examples.models.ex
 | Since :class:`~solver.numerous_solver.NumerousSolver` handles systems of the kind :math:`dy/dt=f(t, y)`, the model
   above has implemented the method `diff` which should be called from :class:`~solver.numerous_solver.NumerousSolver` by
   creating an :class:`~solver.interface.Interface`. Notice how the model is decorated by the
-  :meth:`~solver.interface.Model.with_interface` decorator. If no arguments are provided to the decorator, the interface
-  is assumed to be named ${ModelName}Interface. This is how the model connects to the interface, and it's
-  explained in the next section.
+  :class:`~solver.interface.Model` decorator. Next, we need to create an interface and connect it to this model, which
+   is explained in the next section.
 
 -------------------------------
 Example: Writing an interface
@@ -100,26 +99,42 @@ Example: Writing an interface
 
 If you wish to solve an ODE, you need to implement the :meth:`~solver.interface.Interface.get_deriv` method as well.
 
-| With the example above for :class:`~solver.examples.models.exponential_approach.ExponentialApproach` the correspoding
+| With the example above for :class:`~solver.examples.models.exponential_approach.ExponentialApproach` the corresponding
   interface would be:
 
 .. literalinclude:: ../../src/numerous/solver/examples/models/exponential_approach.py
     :pyobject: ExponentialApproachInterface
     :language: python
 
+Here, we have decorated the interface using the :meth:`~solver.interface.Model.interface` decorator, with arguments
+`model=ExponentialApproach`, where model refers to the
+:class:`~solver.examples.model.exponential_approach.ExponentialApproach` model class. Decorating this way links the
+interface and its model. A model can only have one interface. If we create another interface somewhere in our code, and
+decorate with the same method:
+
+.. Code-Block::
+    @Model.interface(model=ExponentialApproach)
+    class NewInterface(Interface):
+        ...
+
+then the original interface will be overwritten, depending on where in the script/program this is done. A warning will
+be given if a model already has an interface. The warning can be overwritten using the argument `warn=False` in the
+:meth:`~solver.interface.Model.interface` method.
+
 Lets go through the methods of the :class:`~solver.examples.models.exponential_approach.ExponentialApproachInterface`:
 
-* | :meth:`~solver.examples.models.exponential_approach.ExponentialApproachInterface.get_states` returns the states from the
-    model. This is used to save the states in the :class:`~solver.interface.Solution` object.
-* | :meth:`~solver.examples.models.exponential_approach.ExponentialApproachInterface.set_states` takes the states from the
-    solver and transfers them to the model. The solver operates only with :class:`numpy.ndarray` arrays, but the model
-    may use any format to store the states, using dicts, dataframes etc, though jitting always requires that the
+* | :meth:`~solver.examples.models.exponential_approach.ExponentialApproachInterface.get_states` returns the states from
+    the model. This is used to save the states in the :class:`~solver.interface.Solution` object.
+* | :meth:`~solver.examples.models.exponential_approach.ExponentialApproachInterface.set_states` takes the states from
+    the solver and transfers them to the model. The solver operates only with :class:`numpy.ndarray` arrays, but the
+    model may use any format to store the states, using dicts, dataframes etc, though jitting always requires that the
     objects can be understood by `numba <https://numba.pydata.org/numba-doc/latest/index.html>`_. In the example
     the model states are identical to the solver states, so the transfer is straight-forward.
 * | :meth:`~solver.examples.models.exponential_approach.ExponentialApproachInterface.get_deriv` returns the calculated
     derivatives. In the example, the
     :meth:`~solver.examples.models.exponential_approach.ExponentialApproachInterface.get_deriv` calls the model
-    :meth:`~solver.examples.models.exponential_approach.ExponentialApproach.diff` method, and returns the value of this method.
+    :meth:`~solver.examples.models.exponential_approach.ExponentialApproach.diff` method, and returns the value of this
+    method.
 
 | That's it. It may seem more cumbersome than simply specifying the diff function to the solver, such as it is standard
   practise in many ODE solver suites, however the approach allows for much more advanced model to be implemented, and
diff --git a/setup.py b/setup.py
index 9f9d701579d1d2539f57f64618958ed6529d7ffe..b4df1d59ac3e81037a60344272f9af364c74781f 100644
--- a/setup.py
+++ b/setup.py
@@ -28,7 +28,8 @@ setuptools.setup(
     ],
     install_requires=[
         "numba>=0.56.4",
-        "scipy==1.9.1"
+        "scipy==1.9.1",
+        "Deprecated>=1.2.13"
     ],
     packages=setuptools.find_namespace_packages(where='src'),
     package_dir={'': 'src'},
diff --git a/src/numerous/solver/examples/models/bouncing_ball.py b/src/numerous/solver/examples/models/bouncing_ball.py
index 5290f0cb87cec747a5d18d6306f57fb62117335f..e359f72cbf303e3b84587760b5705a4362912737 100644
--- a/src/numerous/solver/examples/models/bouncing_ball.py
+++ b/src/numerous/solver/examples/models/bouncing_ball.py
@@ -131,7 +131,7 @@ import plotly.graph_objects as go
 import plotly
 
 
-@Model.with_interface()
+@Model
 class BouncingBall:
 
     def __init__(self, x0=1, v0=0, f_loss=0.05, g=9.81):
@@ -157,9 +157,9 @@ class BouncingBall:
         return np.array([dxdt, dvdt], dtype='float')
 
 
+@Model.interface(model=BouncingBall)
 class BouncingBallInterface(Interface):
-    def __init__(self, model):
-        self.model: BouncingBall = model
+    model: BouncingBall
 
     def get_deriv(self, t: float, y: np.array) -> np.ascontiguousarray:
         return self.model.diff(t, y)
@@ -185,7 +185,7 @@ class BouncingBallInterface(Interface):
 if __name__ == "__main__":
 
     model = BouncingBall()
-    numsol = NumerousSolver(model=model, use_jit=False)
+    numsol = NumerousSolver(model=model, use_jit=True)
     time = np.append(np.arange(0, 10, 0.1), 10)
     numsol.solve(time)
     t = np.array(numsol.solution.results).T[0, :]
diff --git a/src/numerous/solver/examples/models/controller.py b/src/numerous/solver/examples/models/controller.py
index ed63ce899e1e943fd6b489bceb8802dd2962590e..eeb3c53c00e5dd2df2e551def302394ed4333b90 100644
--- a/src/numerous/solver/examples/models/controller.py
+++ b/src/numerous/solver/examples/models/controller.py
@@ -156,13 +156,13 @@ A note about jitting and the :meth:`Model.component` decorator
 
 The model of the CSTR has been implemented below in the :class:`Tank` model, while the controllers have been implemented
 as :class:`InletFlowController`, and :class:`OutletFlowController`, respectively. The :class:`Tank` model is decorated
-by the :meth:`Model.with_interface` decorator, while the :class:`In/OutletFlowController` models are decorated with the
-:meth:`Model.component` decorator. The reason being there can only be one interface per model, and that is the
-one specified by the :class:`Tank` class, i.e. the :class:`TankInterface`. Since we want to allow the model to be
-jitted by `numba`, we used the :meth:`Model.component` decorator on the :class:`In/OutletFlowController` classes.
-Without this decorator, `numerous solver` cannot be jitted. Allowing jitting is also the reason why the instantiated
-:class:`InletFlowController` and :class:`OutletFlowController` are fed to the :class:`Tank` as an input, since `numba`
-currently cannot infer the type custom objects. Attempting to use jitting would result in an error.
+by the :class:`~solver.interface.Model` decorator, while the :class:`In/OutletFlowController` models are decorated with
+the :meth:`~solver.interface.Model.component` decorator. The reason being there can only be one interface per model,
+and that is the one specified by the :class:`Tank` class, i.e. the :class:`TankInterface`. Since we want to allow the
+model to be jitted by `numba`, we used the :meth:`Model.component` decorator on the :class:`In/OutletFlowController`
+classes. Without this decorator, `numerous solver` cannot be jitted. Allowing jitting is also the reason why the
+instantiated :class:`InletFlowController` and :class:`OutletFlowController` are fed to the :class:`Tank` as an input,
+since `numba` currently cannot infer the type custom objects. Attempting to use jitting would result in an error.
 
 Examples
 --------------
@@ -263,7 +263,7 @@ class InletFlowController:
         self.actual_flow = self.nominal_flow
 
 
-@Model.with_interface()
+@Model
 class Tank:
     def __init__(self, v0: float = 1.0, a: float = 1.0,
                  inlet_flow_controller: InletFlowController = None, outlet_flow_controller: OutletFlowController = None,
@@ -305,11 +305,9 @@ class Tank:
         self.outlet_flow_controller.reset()
         self.initialize()
 
-
+@Model.interface(model=Tank)
 class TankInterface(Interface):
-
-    def __init__(self, model: Tank):
-        self.model = model
+    model: Tank
 
     def set_states(self, y: np.array) -> None:
         y_cstr = y[:4]
diff --git a/src/numerous/solver/examples/models/coupled_tanks.py b/src/numerous/solver/examples/models/coupled_tanks.py
index 0e2a09f26ff54e4f3a33571a4ea3725286ab8cab..03e72b3fb286aba3bcaaeb3796f5e7b75976493b 100644
--- a/src/numerous/solver/examples/models/coupled_tanks.py
+++ b/src/numerous/solver/examples/models/coupled_tanks.py
@@ -66,23 +66,7 @@ from numerous.solver.interface import Interface, Model
 from numerous.solver.numerous_solver import NumerousSolver
 
 
-
-class CoupledTanksInterface(Interface):
-
-    def __init__(self, model):
-        self.model: CoupledTanks = model
-
-    def get_deriv(self, t: float, y: np.array) -> np.array:
-        return self.model.diff(y)
-
-    def set_states(self, y):
-        self.model.y = y
-
-    def get_states(self):
-        return self.model.y
-
-
-@Model.with_interface(interface=CoupledTanksInterface)
+@Model
 class CoupledTanks:
 
     def __init__(self, n_tanks=10, start_volume=10.0, k=1.0):
@@ -120,6 +104,21 @@ class CoupledTanks:
         self.y = np.zeros(self.n_tanks, dtype='float')
         self.y[0] = self.start_volume
 
+
+@Model.interface(model=CoupledTanks)
+class CoupledTanksInterface(Interface):
+    model: CoupledTanks
+
+    def get_deriv(self, t: float, y: np.array) -> np.array:
+        return self.model.diff(y)
+
+    def set_states(self, y):
+        self.model.y = y
+
+    def get_states(self):
+        return self.model.y
+
+
 if __name__ == "__main__":
     n_tanks = 5
     model = CoupledTanks(n_tanks=n_tanks)
diff --git a/src/numerous/solver/examples/models/exponential_approach.py b/src/numerous/solver/examples/models/exponential_approach.py
index 80c1a53a7cfc518a26fda8527fc67ed6f40679ae..26205c3204534b7e272d890eb1c5239454a13877 100644
--- a/src/numerous/solver/examples/models/exponential_approach.py
+++ b/src/numerous/solver/examples/models/exponential_approach.py
@@ -60,7 +60,8 @@ from numerous.solver.numerous_solver import NumerousSolver
 import plotly.graph_objects as go
 import plotly
 
-@Model.with_interface()
+
+@Model
 class ExponentialApproach:
     def __init__(self, k=1.0):
         self.y = np.array([0.0], dtype='float')
@@ -73,10 +74,10 @@ class ExponentialApproach:
         self.y = np.array([0.0], dtype='float')
 
 
+@Model.interface(model=ExponentialApproach)
 class ExponentialApproachInterface(Interface):
 
-    def __init__(self, model: ExponentialApproach):
-        self.model = model
+    model: ExponentialApproach
 
     def get_states(self) -> np.array:
         return self.model.y
diff --git a/src/numerous/solver/examples/models/exponential_approach_no_derivatives.py b/src/numerous/solver/examples/models/exponential_approach_no_derivatives.py
index 1d5b019c4594aab3bc8a10adaffae7e35863d96b..958bd8afcc15c90cd8ead47873795c8a3dc94f3d 100644
--- a/src/numerous/solver/examples/models/exponential_approach_no_derivatives.py
+++ b/src/numerous/solver/examples/models/exponential_approach_no_derivatives.py
@@ -18,7 +18,8 @@ but for now this is not possible with numerous solver.
 Implementation in numerous solver
 ---------------------------------------------------
 
-The :class:`ExponentialApproach` model is wrapped by the :meth:`~solver.interface.Model.with_interface` decorator,
+The :class:`ExponentialApproach` model is wrapped by the :class:`~solver.interface.Model` decorator, while it's
+interface is decorated by the :meth:`~solver.interface.Model.interface` method,
 which links the model to it's interface :class:`ExponentialApproachInterface`. It takes just one parameters as input:
 
 =============   ==============================================================================================
@@ -47,7 +48,6 @@ by calling the :meth:`update_step` and then :meth:`set_states`.
 Examples
 --------------
 Below is an example code that runs the model, and creates a plot for multiple values of :math:`k`.
-
 """
 
 import numpy as np
@@ -59,7 +59,7 @@ import plotly.graph_objects as go
 import plotly
 
 
-@Model.with_interface()
+@Model
 class ExponentialApproach():
     def __init__(self, k=1.0):
         self.y = np.array([0.0])
@@ -72,10 +72,9 @@ class ExponentialApproach():
         self.y = np.array([0.0])
 
 
+@Model.interface(model=ExponentialApproach)
 class ExponentialApproachInterface(Interface):
-
-    def __init__(self, model: ExponentialApproach):
-        self.model = model
+    model: ExponentialApproach
 
     def set_states(self, states):
         self.model.y = states
@@ -88,6 +87,7 @@ class ExponentialApproachInterface(Interface):
         self.set_states(y)
         return SolveEvent.NoneEvent
 
+
 if __name__ == "__main__":
     model = ExponentialApproach(k=1.0)  # Generate the model
     numsol = NumerousSolver(model=model, method='RK45', use_jit=True)  # pass model to numerous solver
diff --git a/src/numerous/solver/interface.py b/src/numerous/solver/interface.py
index 4f9b53ed1ec525fa143a95b613655b9657d1d084..cb954d33e19b2c48f0f814b890e305b41953278e 100644
--- a/src/numerous/solver/interface.py
+++ b/src/numerous/solver/interface.py
@@ -1,5 +1,6 @@
 import inspect
 import logging
+import warnings
 
 import numba as nb
 import numpy as np
@@ -8,6 +9,7 @@ from numba.experimental import jitclass
 from numpy.typing import NDArray
 from typing import Optional
 from abc import ABC
+from deprecated.sphinx import deprecated
 
 from .solve_states import SolveEvent
 
@@ -489,6 +491,8 @@ class Model:
         """
 
         new_class = self.__new__(self.__class__)
+        if hasattr(self.model_class, '_get_interface'):
+            self.interface = self.model_class._get_interface(self.model_class)
         new_class.__init__(self.model_class, self.interface, self.iscomponent)
         new_class._JitHelper = _JitHelper(new_class.model_class, *args, **kwargs)
         new_class.python_model = new_class._JitHelper.python_model
@@ -497,7 +501,30 @@ class Model:
 
         return new_class
 
+
+    @classmethod
+    def interface(cls, model, warn=True):
+        """
+        class method decorator to decorate a model Interface and link the model class to its interface. The model class
+        is monkey patched with the _get_interface() method. Overriding this method can lead to unexpected results.
+        :param model: The model class to link to.
+        :return: wrapped interface as a Model class.
+        """
+        def inner(interface):
+            def get_interface(self):
+                return interface
+            if hasattr(model.model_class, '_get_interface') and warn:
+                logger.warning(f"model {model.model_class} already has a defined interface, which will be overwritten")
+            model.model_class._get_interface = get_interface
+
+            return cls(model.model_class, interface)
+
+        return inner
+
+
     @classmethod
+    @deprecated(version="1.4.0", reason="This method has been deprecated. Please use decorate model with @Model "
+                                        "decorator, and mark interface with @Model.interface(model=model)")
     def with_interface(cls, interface: type(Interface) | str | None = None):
         """decorator to link model to interface. Classes inside the model class must be passed as arguments, if using
         numba compiled solver (jit=True). This is due to limitations of numba.
diff --git a/tests/src/models/exponential_approach_args_and_kwargs.py b/tests/src/models/exponential_approach_args_and_kwargs.py
index 52311a3bb0a5132e6ece78edcbce54b2771265ec..a49f49294c31f0fd8dc8668ce47ac621e61abca0 100644
--- a/tests/src/models/exponential_approach_args_and_kwargs.py
+++ b/tests/src/models/exponential_approach_args_and_kwargs.py
@@ -1,7 +1,7 @@
 import numpy as np
 from numerous.solver.interface import Model, Interface
 
-@Model.with_interface()
+@Model
 class ExponentialApproach:
     def __init__(self, k, y0=0.0):
         self.y = np.array([y0], dtype='float')
@@ -10,12 +10,9 @@ class ExponentialApproach:
     def diff(self, t, y) -> np.array:
         return np.array([self.k * np.exp(-self.k * t)])
 
-
+@Model.interface(model=ExponentialApproach)
 class ExponentialApproachInterface(Interface):
-
-    def __init__(self, model: ExponentialApproach):
-        self.model = model
-
+    model: ExponentialApproach
     def get_states(self) -> np.array:
         return self.model.y
 
diff --git a/tests/src/test_model_and_interface.py b/tests/src/test_model_and_interface.py
index 133cefad9ca8caf73c8a1c2549161cb86c4b4953..c76190ae4fab08903eda8b26345d8df8b8d71990 100644
--- a/tests/src/test_model_and_interface.py
+++ b/tests/src/test_model_and_interface.py
@@ -1,13 +1,17 @@
+import logging
+
 import pytest
 
 import numpy as np
 from conftest import normal_solver, numerous_solver, ABSTOL, RELTOL, get_timerange, RUNMAX
 from numerous.solver.interface import Model as NumerousModel
 from numerous.solver.numerous_solver import NumerousSolver
-from models.exponential_approach_args_and_kwargs import ExponentialApproach as ExponentialApproachArgsAndKwargs
+from models.exponential_approach_args_and_kwargs import ExponentialApproach as ExponentialApproachArgsAndKwargs, \
+    ExponentialApproachInterface
 from models.exponential_approach_numpy_input import ExponentialApproach as ExponentialApproachNumpyInput
 from models.exponential_approach_time_events_list import ExponentialApproach as ExponentialApproachTimeListEvents
 
+
 @pytest.mark.parametrize("solve", [normal_solver])
 @pytest.mark.parametrize("jit", [False, True])
 @pytest.mark.parametrize("Model", [ExponentialApproachArgsAndKwargs])
@@ -55,4 +59,19 @@ def test_init_solver(jit: bool, numerous_solver: numerous_solver, get_timerange:
         solver.reset()
         solver.solve(timerange)
 
-        assert len(solver.solution.time_event_results) == len(time_events[-1]), "expected time events to match solution"
\ No newline at end of file
+        assert len(solver.solution.time_event_results) == len(time_events[-1]), "expected time events to match solution"
+
+@pytest.mark.parametrize("warn", [True, False])
+def test_interface_overwrite(warn: bool, caplog):
+
+    with caplog.at_level(logging.WARNING):
+        model = ExponentialApproachArgsAndKwargs(k=1.0)
+        @NumerousModel.interface(model=ExponentialApproachArgsAndKwargs, warn=warn)
+        class NewInterface(ExponentialApproachInterface):
+            pass
+        if warn:
+            assert "overwritten" in caplog.text, "expected warning"
+        else:
+            assert not "overwritten" in caplog.text, "expected no warning"
+
+
diff --git a/version.txt b/version.txt
index 1892b926767774e9ba91f1e584fa71b4c56abb69..54dd11e46b5877c605270e4b81a70cd4b7fba865 100644
--- a/version.txt
+++ b/version.txt
@@ -1 +1 @@
-1.3.2
+1.4.0-alpha.2