Trying to simulate Weak Value Amplification using QuTiP

Weak Value Amplification is a procedure in which one consider a bipartite Hamiltonian of the form $$H = g A \otimes B$$, where $$A$$ is called system, $$B$$ is called meter and $$g$$ is the extremely small parameter we are interested to measure by amplification. The procedure goes as follows:

(1) The system and meter is initialized by some initial state $$\vert\Psi\rangle=\vert\psi_{i}\rangle \otimes \vert \phi\rangle$$. The evolution of the system $$U = \exp(-igA \otimes B)$$ is approximated to $$1_{A}\otimes 1_{B}-ig A\otimes B$$, since $$g$$ is extremely small. we get

$$U\vert\Psi_{i}\rangle\approx(1_{A}\otimes 1_{B}-ig A\otimes B)\vert\psi_{i}\rangle \otimes \vert \phi\rangle,$$ $$(1_{A}\vert\psi_{i}\rangle\otimes 1_{B}\vert \phi\rangle-ig A\vert\psi_{i}\rangle\otimes B\vert\phi\rangle).$$

(2) After the evolution, the system is post selected onto some final state $$\vert\psi_{f}\rangle$$ (which is nearly orthogonal to the initial state)

$$\langle\psi_{f}\vert U\vert\Psi_{i}\rangle\approx P_{s}(1_{B}\vert\phi\rangle - ig A_{w}B\vert\phi\rangle)$$.

Where $$P_s = \langle\psi_{f}\vert\psi_{i}\rangle$$ is the post-selection probability and $$A_{w} = \langle\psi_{f}\vert A\vert\psi_{i}\rangle/P_{s}$$ is called the weak value.

(3) If i re-exponentiate the above term, we get

$$\frac{\langle\psi_{f}\vert U\vert\Psi_{i}\rangle}{P_{s}}=\exp(-igA_{w}B)\vert\phi\rangle$$

where the exponential on the RHS is the Kraus operator acting on the meter intial state, using which one can estimate the small parameter $$g$$.

My attempt at solving:

(step 1) I prepare the initial state say psi1 and evolved the state according the Hamiltonian using mesolve for a short time t = np.linspace(0, 0.01, 50)

result = mesolve(H, psi1, t, c_ops = [], e_ops =[])

(step 2) The evolved state is result.states, but the problem I am facing is how do I make post-measurement on this with the final state?

One such dumb attempt, I did was taking the final state psif and sandwiched directly with result.states as

eta = psif.dag()*result.states

But my objective was to do the post-selection at a specific time, like say I weakly evolve the system for a short time say 0 to 0.1 seconds and do post-selection exactly at 0.05th second. I am really lost what to do from my (step 1). Is there a way to do post-selection in QuTiP at a specific time?

But QuTiP gives an error about incompatible object sizes if you do this operation, which seems mathematically fine! That is because it prefers that you operate on states in the composite Hilbert space with things defined on the composite Hilbert space. Thus, instead of the shorthand definition of $$\langle\psi_f|U|\Psi_i\rangle$$, we are really using $$\left(\langle\psi_f|\otimes \mathbb{I}\right)U|\Psi_i\rangle.$$ If we define my psif as your psif tensored with the identity on the second Hilbert space, everything will work nicely.
So: the unnormalized post-selected state at time point t[i] can be found through tensor(psif,qeye(N)).dag()*result.states[i], where N is the dimension of the Hilbert space on of the meter. You can verify that this gives psi1 when you $$i=0$$.