diff --git a/doc/changes/2318.feature b/doc/changes/2318.feature
new file mode 100644
index 0000000000..d8b7f54b23
--- /dev/null
+++ b/doc/changes/2318.feature
@@ -0,0 +1,2 @@
+Create `run_from_experiment`, which allows to run stochastic evolution from
+know noise or measurements.
diff --git a/qutip/solver/multitraj.py b/qutip/solver/multitraj.py
index 3d1cf6ee95..032dad3a01 100644
--- a/qutip/solver/multitraj.py
+++ b/qutip/solver/multitraj.py
@@ -234,7 +234,10 @@ def run(self, state, tlist, ntraj=1, *,
     def _initialize_run_one_traj(self, seed, state, tlist, e_ops,
                                  **integrator_kwargs):
         result = self._trajectory_resultclass(e_ops, self.options)
-        generator = self._get_generator(seed)
+        if "generator" in integrator_kwargs:
+            generator = integrator_kwargs.pop("generator")
+        else:
+            generator = self._get_generator(seed)
         self._integrator.set_state(tlist[0], state, generator,
                                    **integrator_kwargs)
         result.add(tlist[0], self._restore_state(state, copy=False))
@@ -287,11 +290,6 @@ def _get_generator(self, seed):
         If the ``seed`` has a ``random`` method, it will be used as the
         generator.
         """
-        if hasattr(seed, 'random'):
-            # We check for the method, not the type to accept pseudo non-random
-            # generator for debug/testing purpose.
-            return seed
-
         if self.options['bitgenerator']:
             bit_gen = getattr(np.random, self.options['bitgenerator'])
             generator = np.random.Generator(bit_gen(seed))
diff --git a/qutip/solver/sode/_noise.py b/qutip/solver/sode/_noise.py
index e1b4592a49..eec07a1485 100644
--- a/qutip/solver/sode/_noise.py
+++ b/qutip/solver/sode/_noise.py
@@ -1,6 +1,6 @@
 import numpy as np
 
-__all__ = ["Wiener"]
+__all__ = ["Wiener", "PreSetWiener"]
 
 
 class Wiener:
@@ -10,31 +10,84 @@ class Wiener:
     def __init__(self, t0, dt, generator, shape):
         self.t0 = t0
         self.dt = dt
-        self.generator = generator
-        self.t_end = t0
         self.shape = shape
-        self.process = np.zeros((1,) + shape, dtype=float)
+        self.generator = generator
+        self.noise = np.zeros((0,) + shape, dtype=float)
+        self.last_W = np.zeros(shape[-1], dtype=float)
+        self.idx_last_0 = 0
 
-    def _extend(self, t):
-        N_new_vals = int((t - self.t_end + self.dt*0.01) // self.dt)
+    def _extend(self, idx):
+        N_new_vals = idx - self.noise.shape[0]
         dW = self.generator.normal(
             0, np.sqrt(self.dt), size=(N_new_vals,) + self.shape
         )
-        W = self.process[-1, :, :] + np.cumsum(dW, axis=0)
-        self.process = np.concatenate((self.process, W), axis=0)
-        self.t_end = self.t0 + (self.process.shape[0] - 1) * self.dt
+        self.noise = np.concatenate((self.noise, dW), axis=0)
 
     def dW(self, t, N):
-        if t + N * self.dt > self.t_end:
-            self._extend(t + N * self.dt)
-        idx0 = int((t - self.t0 + self.dt * 0.01) // self.dt)
-        return np.diff(self.process[idx0:idx0 + N + 1, :, :], axis=0)
+        # Find the index of t.
+        # Rounded to the closest step, but only multiple of dt are expected.
+        idx0 = round((t - self.t0) / self.dt)
+        if idx0 + N - 1 >= self.noise.shape[0]:
+            self._extend(idx0 + N)
+        return self.noise[idx0:idx0 + N, :, :]
 
     def __call__(self, t):
-        if t > self.t_end:
-            self._extend(t)
-        idx = int((t - self.t0 + self.dt * 0.01) // self.dt)
-        return self.process[idx, 0, :]
+        """
+        Return the Wiener process at the closest ``dt`` step to ``t``.
+        """
+        # The Wiener process is not used directly in the evolution, so it's
+        # less optimized than the ``dW`` method.
+
+        # Find the index of t.
+        # Rounded to the closest step, but only multiple of dt are expected.
+        idx = round((t - self.t0) / self.dt)
+        if idx >= self.noise.shape[0]:
+            self._extend(idx + 1)
+
+        if self.idx_last_0 > idx:
+            # Before last call, reseting
+            self.idx_last_0 = 0
+            self.last_W = np.zeros(self.shape[-1], dtype=float)
+
+        self.last_W = self.last_W + np.sum(
+            self.noise[self.idx_last_0:idx+1, 0, :], axis=0
+        )
+
+        self.idx_last_0 = idx
+        return self.last_W
+
+
+class PreSetWiener(Wiener):
+    def __init__(self, noise, tlist, n_sc_ops, heterodyne, is_measurement):
+        if heterodyne:
+            if noise.shape != (n_sc_ops/2, 2, len(tlist)-1):
+                raise ValueError(
+                    "Noise is not of the expected shape: "
+                    f"{(n_sc_ops/2, 2, len(tlist)-1)}"
+                )
+            noise = np.reshape(noise, (n_sc_ops, len(tlist)-1), "C")
+        else:
+            if noise.shape != (n_sc_ops, len(tlist)-1):
+                raise ValueError(
+                    "Noise is not of the expected shape: "
+                    f"{(n_sc_ops, len(tlist)-1)}"
+                )
+
+        self.t0 = tlist[0]
+        self.dt = tlist[1] - tlist[0]
+        self.shape = noise.shape[1:]
+        self.noise = noise.T[:, np.newaxis, :].copy()
+        self.last_W = np.zeros(self.shape[-1], dtype=float)
+        self.idx_last_0 = 0
+        self.is_measurement = is_measurement
+        if self.is_measurement:
+            # Measurements is scaled as <M> + dW / dt
+            self.noise *= self.dt
+            if heterodyne:
+                self.noise /= 2**0.5
+
+    def _extend(self, N):
+        raise ValueError("Requested time is outside the integration range.")
 
 
 class _Noise:
diff --git a/qutip/solver/sode/_sode.pyx b/qutip/solver/sode/_sode.pyx
index d528228d87..a38079686e 100644
--- a/qutip/solver/sode/_sode.pyx
+++ b/qutip/solver/sode/_sode.pyx
@@ -11,9 +11,11 @@ import numpy as np
 
 cdef class Euler:
     cdef _StochasticSystem system
+    cdef bint measurement_noise
 
-    def __init__(self, _StochasticSystem system):
+    def __init__(self, _StochasticSystem system, measurement_noise=False):
         self.system = system
+        self.measurement_noise = measurement_noise
 
     @cython.wraparound(False)
     def run(
@@ -37,9 +39,16 @@ cdef class Euler:
         """
         cdef int i
         cdef _StochasticSystem system = self.system
+        cdef list expect
 
         cdef Data a = system.drift(t, state)
         b = system.diffusion(t, state)
+
+        if self.measurement_noise:
+            expect = system.expect(t, state)
+            for i in range(system.num_collapse):
+                dW[0, i] -= expect[i].real * dt
+
         cdef Data new_state = _data.add(state, a, dt)
         for i in range(system.num_collapse):
             new_state = _data.add(new_state, b[i], dW[0, i])
@@ -72,6 +81,12 @@ cdef class Platen(Euler):
         cdef list d2 = system.diffusion(t, state)
         cdef Data Vt, out
         cdef list Vp, Vm
+        cdef list expect
+
+        if self.measurement_noise:
+            expect = system.expect(t, state)
+            for i in range(system.num_collapse):
+                dW[0, i] -= expect[i].real * dt
 
         out = _data.mul(d1, 0.5)
         Vt = d1.copy()
@@ -106,6 +121,9 @@ cdef class Platen(Euler):
 
 
 cdef class Explicit15(Euler):
+    def __init__(self, _StochasticSystem system):
+        self.system = system
+
     @cython.boundscheck(False)
     @cython.wraparound(False)
     @cython.cdivision(True)
@@ -235,9 +253,11 @@ cdef class Explicit15(Euler):
 
 cdef class Milstein:
     cdef _StochasticSystem system
+    cdef bint measurement_noise
 
-    def __init__(self, _StochasticSystem system):
-        self.system = system
+    def __init__(self, _StochasticSystem system, measurement_noise=False):
+            self.system = system
+            self.measurement_noise = measurement_noise
 
     @cython.wraparound(False)
     def run(self, double t, Data state, double dt, double[:, :, ::1] dW, int ntraj):
@@ -273,6 +293,11 @@ cdef class Milstein:
         iadd_dense(out, state, 1)
         iadd_dense(out, system.a(), dt)
 
+        if self.measurement_noise:
+            expect = system.expect(t, state)
+            for i in range(system.num_collapse):
+                dW[0, i] -= system.expect_i(i).real * dt
+
         for i in range(num_ops):
             iadd_dense(out, system.bi(i), dW[0, i])
 
@@ -289,11 +314,17 @@ cdef class PredCorr:
     cdef Dense euler
     cdef double alpha, eta
     cdef _StochasticSystem system
+    cdef bint measurement_noise
 
-    def __init__(self, _StochasticSystem system, double alpha=0., double eta=0.5):
+    def __init__(
+        self, _StochasticSystem system,
+        double alpha=0., double eta=0.5,
+        measurement_noise=False
+    ):
         self.system = system
         self.alpha = alpha
         self.eta = eta
+        self.measurement_noise = measurement_noise
 
     @cython.wraparound(False)
     def run(self, double t, Data state, double dt, double[:, :, ::1] dW, int ntraj):
@@ -324,6 +355,11 @@ cdef class PredCorr:
 
         system.set_state(t, state)
 
+        if self.measurement_noise:
+            expect = system.expect(t, state)
+            for i in range(system.num_collapse):
+                dW[0, i] -= system.expect_i(i).real * dt
+
         imul_dense(out, 0.)
         iadd_dense(out, state, 1)
         iadd_dense(out, system.a(), dt * (1-alpha))
@@ -350,6 +386,10 @@ cdef class PredCorr:
 
 
 cdef class Taylor15(Milstein):
+    def __init__(self, _StochasticSystem system):
+        self.system = system
+        self.measurement_noise = False
+
     @cython.boundscheck(False)
     @cython.wraparound(False)
     cdef Data step(self, double t, Dense state, double dt, double[:, :] dW, Dense out):
@@ -398,17 +438,17 @@ cdef class Milstein_imp:
     cdef double prev_dt
     cdef dict imp_opt
 
-    def __init__(self, _StochasticSystem system, imp_method=None, imp_options={}):
+    def __init__(self, _StochasticSystem system, solve_method=None, solve_options={}):
         self.system = system
         self.prev_dt = 0
-        if imp_method == "inv":
+        if solve_method == "inv":
             if not self.system.L.isconstant:
                 raise TypeError("The 'inv' integration method requires that the system Hamiltonian or Liouvillian be constant.")
             self.use_inv = True
             self.imp_opt = {}
         else:
             self.use_inv = False
-            self.imp_opt = {"method": imp_method, "options": imp_options}
+            self.imp_opt = {"method": solve_method, "options": solve_options}
 
 
     @cython.wraparound(False)
diff --git a/qutip/solver/sode/itotaylor.py b/qutip/solver/sode/itotaylor.py
index 6ad4df1545..a4541206cd 100644
--- a/qutip/solver/sode/itotaylor.py
+++ b/qutip/solver/sode/itotaylor.py
@@ -17,8 +17,13 @@ class EulerSODE(_Explicit_Simple_Integrator):
 
     - Order: 0.5
     """
+    integrator_options = {
+        "dt": 0.001,
+        "tol": 1e-10,
+    }
     stepper = _sode.Euler
     N_dw = 1
+    _stepper_options = ["measurement_noise"]
 
 
 class Milstein_SODE(_Explicit_Simple_Integrator):
@@ -30,8 +35,13 @@ class Milstein_SODE(_Explicit_Simple_Integrator):
 
     - Order strong 1.0
     """
+    integrator_options = {
+        "dt": 0.001,
+        "tol": 1e-10,
+    }
     stepper = _sode.Milstein
     N_dw = 1
+    _stepper_options = ["measurement_noise"]
 
 
 class Taylor1_5_SODE(_Explicit_Simple_Integrator):
diff --git a/qutip/solver/sode/sode.py b/qutip/solver/sode/sode.py
index e43e8f5ac1..d76b371374 100644
--- a/qutip/solver/sode/sode.py
+++ b/qutip/solver/sode/sode.py
@@ -3,7 +3,7 @@
 from . import _sode
 from ..integrator.integrator import Integrator
 from ..stochastic import StochasticSolver, SMESolver
-from ._noise import Wiener
+from ._noise import Wiener, PreSetWiener
 
 __all__ = ["SIntegrator", "PlatenSODE", "PredCorr_SODE"]
 
@@ -65,7 +65,24 @@ def set_state(self, t, state0, generator):
         """
         self.t = t
         self.state = state0
-        if isinstance(generator, Wiener):
+        stepper_opt = {
+            key: self.options[key]
+            for key in self._stepper_options
+            if key in self.options
+        }
+
+        if isinstance(generator, PreSetWiener):
+            self.wiener = generator
+            if (
+                generator.is_measurement
+                and "measurement_noise" not in self._stepper_options
+            ):
+                raise NotImplementedError(
+                    f"{type(self).__name__} does not support running"
+                    " the evolution from measurements."
+                )
+            stepper_opt["measurement_noise"] = generator.is_measurement
+        elif isinstance(generator, Wiener):
             self.wiener = generator
         else:
             num_collapse = len(self.rhs.sc_ops)
@@ -74,8 +91,8 @@ def set_state(self, t, state0, generator):
                 (self.N_dw, num_collapse)
             )
         self.rhs._register_feedback(self.wiener)
-        opt = [self.options[key] for key in self._stepper_options]
-        self.step_func = self.stepper(self.rhs(self.options), *opt).run
+        rhs = self.rhs(self.options)
+        self.step_func = self.stepper(rhs, **stepper_opt).run
         self._is_set = True
 
     def get_state(self, copy=True):
@@ -228,9 +245,13 @@ class PlatenSODE(_Explicit_Simple_Integrator):
 
     - Order: strong 1, weak 2
     """
-
+    integrator_options = {
+        "dt": 0.001,
+        "tol": 1e-10,
+    }
     stepper = _sode.Platen
     N_dw = 1
+    _stepper_options = ["measurement_noise"]
 
 
 class PredCorr_SODE(_Explicit_Simple_Integrator):
@@ -257,7 +278,7 @@ class PredCorr_SODE(_Explicit_Simple_Integrator):
     }
     stepper = _sode.PredCorr
     N_dw = 1
-    _stepper_options = ["alpha", "eta"]
+    _stepper_options = ["alpha", "eta", "measurement_noise"]
 
     @property
     def options(self):
diff --git a/qutip/solver/sode/ssystem.pxd b/qutip/solver/sode/ssystem.pxd
index a09212c1db..6b78ebe9f8 100644
--- a/qutip/solver/sode/ssystem.pxd
+++ b/qutip/solver/sode/ssystem.pxd
@@ -13,10 +13,13 @@ cdef class _StochasticSystem:
 
     cpdef list diffusion(self, t, Data state)
 
+    cpdef list expect(self, t, Data state)
+
     cpdef void set_state(self, double t, Dense state) except *
 
     cpdef Data a(self)
     cpdef Data bi(self, int i)
+    cpdef complex expect_i(self, int i)
     cpdef Data Libj(self, int i, int j)
     cpdef Data Lia(self, int i)
     cpdef Data L0bi(self, int i)
diff --git a/qutip/solver/sode/ssystem.pyx b/qutip/solver/sode/ssystem.pyx
index be15d4634a..88acefaea4 100644
--- a/qutip/solver/sode/ssystem.pyx
+++ b/qutip/solver/sode/ssystem.pyx
@@ -51,6 +51,12 @@ cdef class _StochasticSystem:
         """
         raise NotImplementedError
 
+    cpdef list expect(self, t, Data state):
+        """
+            Compute the expectation terms for the ``state`` at time ``t``.
+        """
+        raise NotImplementedError
+
     cpdef void set_state(self, double t, Dense state) except *:
         """
             Initialize the set of derrivatives.
@@ -69,6 +75,12 @@ cdef class _StochasticSystem:
         """
         raise NotImplementedError
 
+    cpdef complex expect_i(self, int i):
+        """
+          Expectation value of the ``i``th operator.
+        """
+        raise NotImplementedError
+
     cpdef Data Libj(self, int i, int j):
         """
             bi_n * d bj / dx_n
@@ -144,7 +156,7 @@ cdef class StochasticClosedSystem(_StochasticSystem):
     cpdef list diffusion(self, t, Data state):
         cdef int i
         cdef QobjEvo c_op
-        out = []
+        cdef list out = []
         for i in range(self.num_collapse):
             c_op = self.c_ops[i]
             _out = c_op.matmul_data(t, state)
@@ -153,6 +165,15 @@ cdef class StochasticClosedSystem(_StochasticSystem):
             out.append(_data.add(_out, state, -0.5 * expect))
         return out
 
+    cpdef list expect(self, t, Data state):
+        cdef int i
+        cdef QobjEvo c_op
+        cdef list expect = []
+        for i in range(self.num_collapse):
+            c_op = self.cpcd_ops[i]
+            expect.append(c_op.expect_data(t, state))
+        return expect
+
     def __reduce__(self):
         return (
             StochasticClosedSystem.restore,
@@ -210,7 +231,7 @@ cdef class StochasticOpenSystem(_StochasticSystem):
         cdef int i
         cdef QobjEvo c_op
         cdef complex expect
-        cdef out = []
+        cdef list out = []
         for i in range(self.num_collapse):
             c_op = self.c_ops[i]
             vec = c_op.matmul_data(t, state)
@@ -218,6 +239,16 @@ cdef class StochasticOpenSystem(_StochasticSystem):
             out.append(_data.add(vec, state, -expect))
         return out
 
+    cpdef list expect(self, t, Data state):
+        cdef int i
+        cdef QobjEvo c_op
+        cdef list expect = []
+        for i in range(self.num_collapse):
+            c_op = self.c_ops[i]
+            vec = c_op.matmul_data(t, state)
+            expect.append(_data.trace_oper_ket(vec))
+        return expect
+
     cpdef void set_state(self, double t, Dense state) except *:
         cdef n, l
         self.t = t
@@ -277,6 +308,16 @@ cdef class StochasticOpenSystem(_StochasticSystem):
             self._compute_b()
         return _dense_wrap(self._b[i, :])
 
+    cpdef complex expect_i(self, int i):
+        if not self._is_set:
+            raise RuntimeError(
+                "Derrivatives set for ito taylor expansion need "
+                "to receive the state with `set_state`."
+            )
+        if not self._b_set:
+            self._compute_b()
+        return self.expect_Cv[i]
+
     @cython.boundscheck(False)
     @cython.wraparound(False)
     cdef void _compute_b(self) except *:
@@ -510,6 +551,13 @@ cdef class SimpleStochasticSystem(_StochasticSystem):
             out.append(self.c_ops[i].matmul_data(t, state))
         return out
 
+    cpdef list expect(self, t, Data state):
+        cdef int i
+        cdef list expect = []
+        for i in range(self.num_collapse):
+            expect.append(0j)
+        return expect
+
     cpdef void set_state(self, double t, Dense state) except *:
         self.t = t
         self.state = state
@@ -520,6 +568,9 @@ cdef class SimpleStochasticSystem(_StochasticSystem):
     cpdef Data bi(self, int i):
         return self.c_ops[i].matmul_data(self.t, self.state)
 
+    cpdef complex expect_i(self, int i):
+        return 0j
+
     cpdef Data Libj(self, int i, int j):
         bj = self.c_ops[i].matmul_data(self.t, self.state)
         return self.c_ops[j].matmul_data(self.t, bj)
diff --git a/qutip/solver/stochastic.py b/qutip/solver/stochastic.py
index ee35c63f57..9a81fcf775 100644
--- a/qutip/solver/stochastic.py
+++ b/qutip/solver/stochastic.py
@@ -1,6 +1,7 @@
 __all__ = ["smesolve", "SMESolver", "ssesolve", "SSESolver"]
 
 from .sode.ssystem import StochasticOpenSystem, StochasticClosedSystem
+from .sode._noise import PreSetWiener
 from .result import MultiTrajResult, Result, ExpectOp
 from .multitraj import _MultiTrajRHS, MultiTrajSolver
 from .. import Qobj, QobjEvo
@@ -9,28 +10,28 @@
 from functools import partial
 from .solver_base import _solver_deprecation
 from ._feedback import _QobjFeedback, _DataFeedback, _WienerFeedback
+from time import time
 
 
 class StochasticTrajResult(Result):
     def _post_init(self, m_ops=(), dw_factor=(), heterodyne=False):
         super()._post_init()
-        self.W = []
-        self.m_ops = []
-        self.m_expect = []
-        self.dW_factor = dw_factor
+        self.noise = []
         self.heterodyne = heterodyne
-        for op in m_ops:
-            f = self._e_op_func(op)
-            self.W.append([0.0])
-            self.m_expect.append([])
-            self.m_ops.append(ExpectOp(op, f, self.m_expect[-1].append))
-            self.add_processor(self.m_ops[-1]._store)
+        if self.options["store_measurement"]:
+            self.m_ops = []
+            self.m_expect = []
+            self.dW_factor = dw_factor
+            for op in m_ops:
+                f = self._e_op_func(op)
+                self.m_expect.append([])
+                self.m_ops.append(ExpectOp(op, f, self.m_expect[-1].append))
+                self.add_processor(self.m_ops[-1]._store)
 
     def add(self, t, state, noise=None):
         super().add(t, state)
-        if noise is not None and self.options["store_measurement"]:
-            for i, dW in enumerate(noise):
-                self.W[i].append(self.W[i][-1] + dW)
+        if noise is not None:
+            self.noise.append(noise)
 
     @property
     def wiener_process(self):
@@ -43,7 +44,11 @@ def wiener_process(self):
             (len(sc_ops), 2, len(tlist))
         for heterodyne detection.
         """
-        W = np.array(self.W)
+        W = np.zeros(
+            (self.noise[0].shape[0], len(self.times)),
+            dtype=np.float64
+        )
+        np.cumsum(np.array(self.noise).T, axis=1, out=W[:, 1:])
         if self.heterodyne:
             W = W.reshape(-1, 2, W.shape[1])
         return W
@@ -59,10 +64,10 @@ def dW(self):
             (len(sc_ops), 2, len(tlist)-1)
         for heterodyne detection.
         """
-        dw = np.diff(self.W, axis=1)
+        noise = np.array(self.noise).T
         if self.heterodyne:
-            dw = dw.reshape(-1, 2, dw.shape[1])
-        return dw
+            return noise.reshape(-1, 2, noise.shape[1])
+        return noise
 
     @property
     def measurement(self):
@@ -75,15 +80,30 @@ def measurement(self):
             (len(sc_ops), 2, len(tlist)-1)
         for heterodyne detection.
         """
-        dts = np.diff(self.times)
-        m_expect = np.array(self.m_expect)[:, 1:]
-        noise = np.einsum(
-            "i,ij,j->ij", self.dW_factor, np.diff(self.W, axis=1), (1 / dts)
+        if not self.options["store_measurement"]:
+            return None
+        elif self.options["store_measurement"] == "start":
+            m_expect = np.array(self.m_expect)[:, :-1]
+        elif self.options["store_measurement"] == "middle":
+            m_expect = np.apply_along_axis(
+                lambda m: np.convolve(m, [0.5, 0.5], "valid"),
+                axis=1, arr=self.m_expect,
+            )
+        elif self.options["store_measurement"] in ["end", True]:
+            m_expect = np.array(self.m_expect)[:, 1:]
+        else:
+            raise ValueError(
+                "store_measurement must be in {'start', 'middle', 'end', ''}, "
+                f"not {self.options['store_measurement']}"
+            )
+        noise = np.array(self.noise).T
+        noise_scaled = np.einsum(
+            "i,ij,j->ij", self.dW_factor, noise, (1 / np.diff(self.times))
         )
         if self.heterodyne:
             m_expect = m_expect.reshape(-1, 2, m_expect.shape[1])
-            noise = noise.reshape(-1, 2, noise.shape[1])
-        return m_expect + noise
+            noise_scaled = noise_scaled.reshape(-1, 2, noise_scaled.shape[1])
+        return m_expect + noise_scaled
 
 
 class StochasticResult(MultiTrajResult):
@@ -319,9 +339,10 @@ def smesolve(
           | Whether or not to store the state vectors or density matrices.
             On `None` the states will be saved if no expectation operators are
             given.
-        - | store_measurement: bool
-          | Whether to store the measurement and wiener process for each
-            trajectories.
+        - | store_measurement: str, {'start', 'middle', 'end', ''}
+          | Whether and how to store the measurement for each trajectories.
+            'start', 'middle', 'end' indicate when in the interval the
+            expectation value of the ``m_ops`` is taken.
         - | keep_runs_results : bool
           | Whether to store results from all trajectories or just store the
             averages.
@@ -442,9 +463,10 @@ def ssesolve(
           | Whether or not to store the state vectors or density matrices.
             On `None` the states will be saved if no expectation operators are
             given.
-        - | store_measurement: bool
-            Whether to store the measurement and wiener process, or brownian
-            noise for each trajectories.
+        - | store_measurement: str, {'start', 'middle', 'end', ''}
+          | Whether and how to store the measurement for each trajectories.
+            'start', 'middle', 'end' indicate when in the interval the
+            expectation value of the ``m_ops`` is taken.
         - | keep_runs_results : bool
           | Whether to store results from all trajectories or just store the
             averages.
@@ -514,7 +536,7 @@ class StochasticSolver(MultiTrajSolver):
         "num_cpus": None,
         "bitgenerator": None,
         "method": "platen",
-        "store_measurement": False,
+        "store_measurement": "",
     }
 
     def _trajectory_resultclass(self, e_ops, options):
@@ -626,6 +648,97 @@ def _integrate_one_traj(self, seed, tlist, result):
             result.add(t, self._restore_state(state, copy=False), noise)
         return seed, result
 
+    def run_from_experiment(
+        self, state, tlist, noise, *,
+        args=None, e_ops=(), measurement=False,
+    ):
+        """
+        Run a single trajectory from a given state and noise.
+
+        Parameters
+        ----------
+        state : Qobj
+            Initial state of the system.
+
+        tlist : array_like
+            List of times for which to evaluate the state. The tlist must
+            increase uniformly.
+
+        noise : array_like
+            Noise for each time step and each stochastic collapse operators.
+            For homodyne detection, ``noise[i, t_idx]`` is the Wiener
+            increments between ``tlist[t_idx]`` and ``tlist[t_idx+1]`` for the
+            i-th sc_ops.
+            For heterodyne detection, an extra dimension is added for the pair
+            of measurement: ``noise[i, j, t_idx]``with ``j`` in ``{0,1}``.
+
+        args : dict, optional
+            Arguments to pass to the Hamiltonian and collapse operators.
+
+        e_ops : list, optional
+            List of operators for which to evaluate expectation values.
+
+        measurement : bool, default : False
+            Whether the passed noise is the Wiener increments ``dW`` (gaussian
+            noise with standard derivation of dt**0.5), or the measurement.
+
+            Homodyne measurement is::
+
+              noise[i][t] = dW/dt + expect(sc_ops[i] + sc_ops[i].dag, state[t])
+
+            Heterodyne measurement is::
+
+              noise[i][0][t] = dW/dt * 2**0.5
+                + expect(sc_ops[i] + sc_ops[i].dag, state[t])
+
+              noise[i][1][t] = dW/dt * 2**0.5
+                -1j * expect(sc_ops[i] - sc_ops[i].dag, state[t])
+
+            Note that this function expects the expectation values to be taken
+            at the start of the time step, corresponding to the "start" setting
+            for the "store_measurements" option.
+
+            Only available for limited integration methods.
+
+        Returns
+        -------
+        result : StochasticTrajResult
+            Result of the trajectory.
+
+        Notes
+        -----
+        Only default values of `m_ops` and `dW_factors` are supported.
+        """
+        start_time = time()
+        self._argument(args)
+        stats = self._initialize_stats()
+        dt = tlist[1] - tlist[0]
+        if not np.allclose(dt, np.diff(tlist)):
+            raise ValueError("tlist must be evenly spaced.")
+        generator = PreSetWiener(
+            noise, tlist, len(self.rhs.sc_ops), self.heterodyne, measurement
+        )
+        state0 = self._prepare_state(state)
+        try:
+            old_dt = None
+            if "dt" in self._integrator.options:
+                old_dt = self._integrator.options["dt"]
+                self._integrator.options["dt"] = dt
+            mid_time = time()
+            result = self._initialize_run_one_traj(
+                None, state0, tlist, e_ops, generator=generator
+            )
+            _, result = self._integrate_one_traj(None, tlist, result)
+        except Exception as err:
+            if old_dt is not None:
+                self._integrator.options["dt"] = old_dt
+            raise
+
+        stats['preparation time'] += mid_time - start_time
+        stats['run time'] = time() - mid_time
+        result.stats.update(stats)
+        return result
+
     @classmethod
     def avail_integrators(cls):
         if cls is StochasticSolver:
@@ -649,8 +762,10 @@ def options(self):
             On `None` the states will be saved if no expectation operators are
             given.
 
-        store_measurement: bool, default: False
-            Whether to store the measurement for each trajectories.
+        store_measurement: str, {'start', 'middle', 'end', ''}, default: ""
+            Whether and how to store the measurement for each trajectories.
+            'start', 'middle', 'end' indicate when in the interval the
+            expectation value of the ``m_ops`` is taken.
             Storing measurements will also store the wiener process, or
             brownian noise for each trajectories.
 
@@ -795,7 +910,7 @@ class SMESolver(StochasticSolver):
         "num_cpus": None,
         "bitgenerator": None,
         "method": "platen",
-        "store_measurement": False,
+        "store_measurement": "",
     }
 
 
@@ -839,5 +954,5 @@ class SSESolver(StochasticSolver):
         "num_cpus": None,
         "bitgenerator": None,
         "method": "platen",
-        "store_measurement": False,
+        "store_measurement": "",
     }
diff --git a/qutip/tests/solver/test_sode_method.py b/qutip/tests/solver/test_sode_method.py
index 356a2ce70a..5416f5335a 100644
--- a/qutip/tests/solver/test_sode_method.py
+++ b/qutip/tests/solver/test_sode_method.py
@@ -60,7 +60,7 @@ def _make_oper(kind, N):
     pytest.param("Euler", 0.5, {}, id="Euler"),
     pytest.param("Milstein", 1.0, {}, id="Milstein"),
     pytest.param("Milstein_imp", 1.0, {}, id="Milstein implicit"),
-    pytest.param("Milstein_imp", 1.0, {"imp_method": "inv"},
+    pytest.param("Milstein_imp", 1.0, {"solve_method": "inv"},
                  id="Milstein implicit inv"),
     pytest.param("Platen", 1.0, {}, id="Platen"),
     pytest.param("PredCorr", 1.0, {}, id="PredCorr"),
@@ -68,7 +68,7 @@ def _make_oper(kind, N):
     pytest.param("Taylor15", 1.5, {}, id="Taylor15"),
     pytest.param("Explicit15", 1.5, {}, id="Explicit15"),
     pytest.param("Taylor15_imp", 1.5, {}, id="Taylor15 implicit"),
-    pytest.param("Taylor15_imp", 1.5, {"imp_method": "inv"},
+    pytest.param("Taylor15_imp", 1.5, {"solve_method": "inv"},
                  id="Taylor15 implicit inv"),
 ])
 @pytest.mark.parametrize(['H', 'sc_ops'], [
@@ -79,7 +79,7 @@ def _make_oper(kind, N):
     pytest.param("qeye", ["qeye", "destroy", "destroy2"], id='3 sc_ops'),
 ])
 def test_methods(H, sc_ops, method, order, kw):
-    if kw == {"imp_method": "inv"} and ("td" in H or "td" in sc_ops[0]):
+    if kw == {"solve_method": "inv"} and ("td" in H or "td" in sc_ops[0]):
         pytest.skip("inverse method only available for constant cases.")
     N = 5
     H = _make_oper(H, N)
diff --git a/qutip/tests/solver/test_stochastic.py b/qutip/tests/solver/test_stochastic.py
index b692312c38..eaf2f9c015 100644
--- a/qutip/tests/solver/test_stochastic.py
+++ b/qutip/tests/solver/test_stochastic.py
@@ -281,10 +281,43 @@ def test_reuse_seeds():
 
 
 @pytest.mark.parametrize("heterodyne", [True, False])
-def test_m_ops(heterodyne):
+def test_measurements(heterodyne):
     N = 10
     ntraj = 1
 
+    H = num(N)
+    sc_ops = [destroy(N)]
+    psi0 = basis(N, N-1)
+
+    times = np.linspace(0, 1.0, 11)
+
+    solver = SMESolver(H, sc_ops, heterodyne=heterodyne)
+
+    solver.options["store_measurement"] = "start"
+    res_start = solver.run(psi0, times, ntraj=ntraj, seeds=1)
+
+    solver.options["store_measurement"] = "middle"
+    res_middle = solver.run(psi0, times, ntraj=ntraj, seeds=1)
+
+    solver.options["store_measurement"] = "end"
+    res_end = solver.run(psi0, times, ntraj=ntraj, seeds=1)
+
+    diff = np.sum(np.abs(res_end.measurement[0] - res_start.measurement[0]))
+    assert diff > 0.1 # Each measurement should be different by ~dt
+    np.testing.assert_allclose(
+        res_middle.measurement[0] * 2,
+        res_start.measurement[0] + res_end.measurement[0],
+    )
+
+    np.testing.assert_allclose(
+        np.diff(res_start.wiener_process[0][0]), res_start.dW[0][0]
+    )
+
+
+@pytest.mark.parametrize("heterodyne", [True, False])
+def test_m_ops(heterodyne):
+    N = 10
+
     H = num(N)
     sc_ops = [destroy(N), qeye(N)]
     psi0 = basis(N, N-1)
@@ -294,7 +327,7 @@ def test_m_ops(heterodyne):
 
     times = np.linspace(0, 1.0, 51)
 
-    options = {"store_measurement": True,}
+    options = {"store_measurement": "end",}
 
     solver = SMESolver(H, sc_ops, heterodyne=heterodyne, options=options)
     solver.m_ops = m_ops
@@ -322,7 +355,6 @@ def test_m_ops(heterodyne):
 
 
 def test_feedback():
-    tol = 0.05
     N = 10
     ntraj = 2
 
@@ -345,7 +377,9 @@ def func(t, A, W):
     solver = SMESolver(H, sc_ops=sc_ops, heterodyne=False, options=options)
     results = solver.run(psi0, times, e_ops=[num(N)], ntraj=ntraj)
 
-    assert np.all(results.expect[0] > 6.-1e-6)
+    # If this was deterministic, it should never go under `6`.
+    # We add a tolerance ~dt due to the stochatic part.
+    assert np.all(results.expect[0] > 6. - 0.001)
     assert np.all(results.expect[0][-20:] < 6.7)
 
 
@@ -380,3 +414,91 @@ def test_small_step_warnings(method):
             qeye(2), basis(2), [0, 0.0000001], [qeye(2)],
             options={"method": method}
         )
+
+
+@pytest.mark.parametrize("method", ["euler", "platen"])
+@pytest.mark.parametrize("heterodyne", [True, False])
+def test_run_from_experiment_close(method, heterodyne):
+    N = 10
+
+    H = num(N)
+    a = destroy(N)
+    sc_ops = [a, a @ a + (a @ a).dag()]
+    psi0 = basis(N, N-1)
+    tlist = np.linspace(0, 0.1, 251)
+    options = {
+        "store_measurement": "start",
+        "dt": tlist[1],
+        "store_states": True,
+        "method": method,
+    }
+    solver = SSESolver(H, sc_ops, heterodyne, options=options)
+    res_forward = solver.run(psi0, tlist, 1, e_ops=[H])
+    res_backward = solver.run_from_experiment(
+        psi0, tlist, res_forward.dW[0], e_ops=[H]
+    )
+    res_measure = solver.run_from_experiment(
+        psi0, tlist, res_forward.measurement[0], e_ops=[H], measurement=True
+    )
+
+    np.testing.assert_allclose(
+        res_backward.measurement, res_forward.measurement[0], atol=1e-10
+    )
+    np.testing.assert_allclose(
+        res_measure.measurement, res_forward.measurement[0], atol=1e-10
+    )
+
+    np.testing.assert_allclose(res_backward.dW, res_forward.dW[0], atol=1e-10)
+    np.testing.assert_allclose(res_measure.dW, res_forward.dW[0], atol=1e-10)
+
+    np.testing.assert_allclose(
+        res_backward.expect, res_forward.expect, atol=1e-10
+    )
+    np.testing.assert_allclose(
+        res_measure.expect, res_forward.expect, atol=1e-10
+    )
+
+
+@pytest.mark.parametrize(
+    "method", ["euler", "milstein", "platen", "pred_corr"]
+)
+@pytest.mark.parametrize("heterodyne", [True, False])
+def test_run_from_experiment_open(method, heterodyne):
+    N = 10
+
+    H = num(N)
+    a = destroy(N)
+    sc_ops = [a, a.dag() * 0.1]
+    psi0 = basis(N, N-1)
+    tlist = np.linspace(0, 1, 251)
+    options = {
+        "store_measurement": "start",
+        "dt": tlist[1],
+        "store_states": True,
+        "method": method,
+    }
+    solver = SMESolver(H, sc_ops, heterodyne, options=options)
+    res_forward = solver.run(psi0, tlist, 1, e_ops=[H])
+    res_backward = solver.run_from_experiment(
+        psi0, tlist, res_forward.dW[0], e_ops=[H]
+    )
+    res_measure = solver.run_from_experiment(
+        psi0, tlist, res_forward.measurement[0], e_ops=[H], measurement=True
+    )
+
+    np.testing.assert_allclose(
+        res_backward.measurement, res_forward.measurement[0], atol=1e-10
+    )
+    np.testing.assert_allclose(
+        res_measure.measurement, res_forward.measurement[0], atol=1e-10
+    )
+
+    np.testing.assert_allclose(res_backward.dW, res_forward.dW[0], atol=1e-10)
+    np.testing.assert_allclose(res_measure.dW, res_forward.dW[0], atol=1e-10)
+
+    np.testing.assert_allclose(
+        res_backward.expect, res_forward.expect, atol=1e-10
+    )
+    np.testing.assert_allclose(
+        res_measure.expect, res_forward.expect, atol=1e-10
+    )