So, I am about 6 years late to the punch, but here is a computational experiment setup. I'd like to code this up some day when I have enough time, it's quite doable.
Theory
You have two Hilbert spaces, one whose elements are the usual wavefunctions $\psi(x,y)\in\mathcal{H}_x$, and one which is your detector which says whether a particle passed through a left or right slit. Say your detector starts in state $(|L\rangle+|R\rangle)/\sqrt{2}$: a 50% superposition of the left and right slits. I'll denote this Hilbert space $\mathcal{H}_d=\operatorname{span}(|L\rangle,|R\rangle)$
The tensor product space of the particle and detector will be described by a function $\psi(x,y,d)$ where $d$ is one of $L$ or $R$. So for example if we had our initial wavepacket of our electron as $\psi_0(x,y)$ and the 50/50 detector state above, our initial full wavefunction would be $\psi(x,y,L)=\frac{1}{\sqrt{2}}\psi_0(x,y)$ and $\psi(x,y,R)=\frac{1}{\sqrt{2}}\psi_0(x,y)$.
For the Hamiltonian, the potential barrier acts as a potential on $\mathcal{H}_x$ and as identity on $\mathcal{H}_d$. So we have a $V(x,y)$ which is infinite except for two blocks cut out: the slits. I'll write the identity on $\mathcal{H}_d$ as $\mathbf{1}_d$.
To model the detection of particles, you can have a test function $h(x,y)$ which is $1$ inside the left slit, $-1$ inside the right slit, and $0$ elsewhere. If a particle is detected, we want to move the amplitude for detection to the correct state $|L\rangle$ or $|R\rangle$, so we should add a term something like the Pauli $\sigma_y$ matrix*. The full Hamiltonian will look like:
$$H=(p_x^2+p_y^2+V(x,y))\otimes \mathbf{1}_d + I h(x,y)\otimes \sigma_y$$
Where $I$ is the interaction strength.
Simulation
For numerical simulation, this turns into the following PDE ($m=\hbar=1$ WLOG):
\begin{align*}
i \partial_t \psi(x,y,L,t)&=(-\frac{1}{2}\partial_x^2-\frac{1}{2}\partial_y^2+V(x,y))\psi(x,y,L,t) -i I h(x,y)\psi(x,y,R,t)\\
i \partial_t \psi(x,y,R,t)&=(-\frac{1}{2}\partial_x^2-\frac{1}{2}\partial_y^2+V(x,y))\psi(x,y,R,t) +i I h(x,y)\psi(x,y,L,t)
\end{align*}
You are basically simulating two separate wavefunctions, one for $L$ and one for $R$, but the term $I$ can mix the two wavefunctions and will do so in a way that disturbs interference. If $\psi(x,y,L)$ is plotted in blue and $\psi(x,y,R)$ is plotted in red, then for $I=0$ we see a purple pattern showing interference. For $I$ tuned to a nice value, we'll see the interference pattern disappearing, and a blue-purple-red bump with no interference pattern showing.
Comment on Interactions
The usual objection (also mentioned in Feynman) is that "well we're turning on an interaction $I$, it's no surprise that messes with the wavefunction." The point is that there's no way whatsoever to get position information out of the Hilbert space of the particle $\mathcal{H}_x$ and into the Hilbert space of the detector $\mathcal{H}_d$ without destroying the interference pattern. None, nada, zilch.
*I choose $\sigma_y$ because I know this will cause states to evolve like $e^{-i H t}$, and $-i \sigma_y t$ is $t\begin{bmatrix} 0 & -1 \\ 1 & 0 \end{bmatrix}$, which is a rotation matrix and is easier for me to visualize.
