If trying to use the wave function description, one cannot always uniquely determine the wave function of the sub-system qubit 1. This roots from the famous EPR (Einstein-Podolsky-Rosen, 1935) paradox.
To understand this, considering the following entangled state between qubit 1 and qubits 2-3
$$
\begin{align}
|\psi \rangle_{1, 23}
& = \frac{1}{\sqrt{2}}(|0\rangle|00\rangle + |1\rangle|11\rangle) \\
& = \frac{1}{\sqrt{2}}(|+\rangle|++\rangle + |-\rangle|--\rangle)
\end{align}
$$
where
$$
|+\rangle \equiv \frac{|0\rangle + |1\rangle}{\sqrt{2}} \\
|-\rangle \equiv \frac{|0\rangle - |1\rangle}{\sqrt{2}} \\
|++\rangle \equiv \frac{|00\rangle + |11\rangle}{\sqrt{2}} \\
|--\rangle \equiv \frac{|00\rangle - |11\rangle}{\sqrt{2}}
$$
If one chooses to measure the qubits 2-3 in {$|00\rangle, |11\rangle$}, the qubit 1 would be produced in {$|0\rangle, |1\rangle$}. Alternatively, if one chooses to measure the qubits 2-3 in {$|++\rangle, |--\rangle$}, the qubit 1 would be produced in {$|+\rangle, |-\rangle$}. But, {$|0\rangle, |1\rangle$} and {$|+\rangle, |-\rangle$} are two different wave functional basis. This paradox shows that the wave function for the sub-system qubit 1 cannot be uniquely written out. It is dependent on the measurement choice to qubits 2-3.
The proper way to express the state of the subsystem qubit 1 is to adopt the density operator formulation. Essentially the density operator was developed to describe the mixed, uncertain, quantum wave functions. For the pure 3-qubit quantum state, it is simply written as
$$ \hat{\rho}_{1, 2, 3} \equiv |\psi(1, 2, 3)\rangle \langle \psi(1, 2, 3)|$$
The uncertain state for the sub-system qubit 1 is obtained by tracing over the sub-system basis of qubits 2 and 3.
$$
\hat{\rho}_{1}
= Tr_{2,3} \left[ |\psi(1, 2, 3)\rangle \langle \psi(1, 2, 3)| \right]
$$
More explicitly
$$
\hat{\rho}_{1} = \sum_{i_2, i_3 \in \{0, 1\}}
{\langle i_2, i_3| \left[ |\psi(1, 2, 3)\rangle \langle \psi(1, 2, 3)| \right] |i_2, i_3\rangle }
\space \space (1)
$$
Applying equation (1) to your example, $ U|\psi\rangle \equiv (a \vert{0}\rangle +b \vert{1}\rangle)\otimes\vert{01}\rangle$, the density operator for qubit 1 would be
$$
\hat{\rho}_{1}(U|\psi\rangle) = (b |1\rangle +a |0\rangle)(a \langle 0| +b \langle 1|)
$$
So your example simply demonstrates a pure wave function state for qubit 1, because qubits 2-3 isn't really entangled with qubit 1.
Similarly, applying equation (1) to the entangled example mentioned above, $|\psi \rangle_{1, 23} \equiv \frac{1}{\sqrt{2}}(|0\rangle|00\rangle + |1\rangle|11\rangle)$, the density operator for qubit 1 would be
$$
\begin{align}
\hat{\rho}_{1}(|\psi \rangle_{1, 23}) &= \frac{1}{2}|0\rangle \langle 0| + \frac{1}{2}|1\rangle \langle 1| \\
&= \frac{1}{2}|+\rangle \langle +| + \frac{1}{2}|-\rangle \langle -|
\end{align}
$$
This example demonstrates a mixed, uncertain quantum state of qubit 1, as described above.
[Answer to EDIT]
As to the recipe, or an algorithm, to resolve the qubit-1 state and code the comparison function, I came up with a heuristic one. The idea is simple.
The disentangled state for qubit-1 and qubit-2-3 can be written as
$$
\begin{align}
|\psi_1\rangle |\psi_{23}\rangle &\equiv (a'|0\rangle + b'|1\rangle)
(c|00\rangle +d|01\rangle +e|10\rangle +f|11\rangle) \\
&= a'|0\rangle (c|00\rangle +d|01\rangle +e|10\rangle +f|11\rangle) \\
&+ b'|1\rangle (c|00\rangle +d|01\rangle +e|10\rangle +f|11\rangle)
\end{align}
$$
which is simply spanned over the 8-basis of the 3 qubits. The first 4 terms and the last 4 terms contains the reduced two wave functions of qubit-2-3. The coefficients of these vectors are linearly dependent (i.e., another way to say they are the same wave functions). So we can encode any 3-qubit state in an 8-complex-component vector
$$
[c_0, c_1, c_2, c_3, \space\space c_4, c_5, c_6, c_7]
$$
The state is disentangled iff the first 4 terms has some common factor ($\frac{a'}{b'}$) to the last 4 terms
$$
\frac{a'}{b'} \equiv c' = \frac{c_0}{c_4} = \frac{c_1}{c_5} = \frac{c_2}{c_6} = \frac{c_3}{c_7}
$$
Just make sure the first and the last 4 terms of the encoded vector correspond to qubit-1's $|0\rangle$ and $|1\rangle$ states. In your python code, you need to handle the edge cases when any of denominators equal to zero.
If step 1 returns disentanglement true and the common complex factor c', the qubit-1's state is immediately retrieved by
$$
|\psi_1\rangle \propto c'|0\rangle + |1\rangle
$$
Therefore the comparison function for qubit 1's state returns true if
$$
c' == \frac{a}{b}
$$
assuming the original copy of qubit 1 is given by $a|0\rangle + b|1\rangle$. Again, you need to handle the case when the denominator equals to zero in the python code.