Extract single qubit state from combined state in QuTiP using a Hamiltonian This questions is very similar to another question asked a couple of months ago. I would like to extract the state of single qubits from combined states, using QuTip (the Quantum Toolbox in Python) using a Hamiltonian simulator.
Example
For example, suppose I begin with two qubits, |0⟩ and |1⟩. I can create these in QuTiP (as q1 and q2) as follows:
q1 = basis(2,1)
q2 = basis(2,0)

$q_1 = \begin{bmatrix}0\\1\end{bmatrix} \quad q_2 = \begin{bmatrix}1\\0\end{bmatrix}$
I can then combine q1 and q2 using tensor multiplication:
q1q2 = tensor(q1, q2)

$q_1 \otimes q_2 = \begin{bmatrix}0\\0\\1\\0\end{bmatrix}$
Now, suppose that I apply a CNOT gate to $q_1\otimes q_2$. We can define CNOT as:
circ =  QubitCircuit(2)
circ.add_gate("CNOT",targets=[1],controls=[0])
U = gate_sequence_product(circ.propagators())

$U=\begin{bmatrix}
1&0&0&0 \\
0&1&0&0 \\
0&0&0&1 \\
0&0&1&0 \\
\end{bmatrix}$
Probability of $q_2$
I now want to figure what q2 would look like after application of this gate. We can do this using the partial trace of the density matrix:
M = ket2dm(U*q1q2).ptrace(0)
q_2' = np.diag(np.sqrt(M.full()))

$M = \begin{bmatrix}0&0\\0&1 \end{bmatrix}\\
q_2' = \begin{bmatrix}0\\1 \end{bmatrix}$
Probability of $q_2$ version 2
I would like to do the same now, but using a Hamiltonian simulator. For this, I use the built-in master equation solver which solves:
$iℏ\frac{d}{d}t|ψ⟩ = H |ψ⟩$
The code for this is:
times = np.linspace(0.0, 20.0, 100)
data = mesolve(U, q1q2, times)
M2 = ket2dm(data.states[-1]).ptrace(0)
q_2'' = np.diag(np.sqrt(M2.full()))

$M_2 = \begin{bmatrix}0.167&0.373j\\-0.373j&0.833 \end{bmatrix}\\
q_2'' = \begin{bmatrix}0.40809641\\0.91293884 \end{bmatrix}$
Which is not what I expected. Question: can I simply use my unitary matrix here or do I need to do some transformation before/after running the simulator? Are there other ways to come up with the probabilities of one qubit, besides using ptrace?
 A: Say $U$ is a unitary evolution (the CNOT in your case), and $\vert\psi\rangle$ an initial state.
Evolving $|\psi\rangle$ with $U$ amounts to computing $U\lvert\psi\rangle$.
If on the other hand you want to study the dynamics of this evolution, then you need a generator for the evolution. In other words, you want an Hamiltonian $H$ such that after a certain time $t$ you get $$e^{-i t H}\lvert\psi\rangle=U\lvert\psi\rangle.$$
Such an Hamiltonian is what you need to feed to mesolve. Here is a working version of the code doing just this:
import numpy as np
import scipy
import qutip

cnot_generator = qutip.Qobj(1j * scipy.linalg.logm(qutip.cnot().full()), dims=[[2] * 2] * 2)
initial_state = qutip.tensor(qutip.basis(2, 1), qutip.basis(2, 0))
times = np.arange(0, 1.1, 0.1)
qutip.mesolve(cnot_generator, initial_state, times).states

Note how the last state given by the above is (in good approximation) the same you obtained yourself.
I used ket notation here because, being everything unitary, there is no need to use density matrices. You may adapt the code appropriately to use density matrices.
