For 1: It does not get just one final state, it gets a mapping from all basis states to all corresponding final states, which is one way to define a superoperator (just like how a unitary operator is equivalently defined by how it maps all basis states).
For 3: The Clifford approach might not work, but I suspect the tensor network approach will still be pretty fruitful.
Concerning your question about setting it up as an ODE: Yes, internally everything done by QuantumOptics.jl is to convert the nice complicated unitaries or states (with all their basis and tensor product structure) into simple matrices and feed these to solvers from DifferentialEquations.jl. If you want to, you can use QuantumOptics.jl as a convenient way to prepare these matrices and then feed everything to DifferentialEquations.jl yourself.
E.g. after you have created whatever QuantumOptics.jl you like, you can extract the underlying matrix from the data attribute. Here is an example where the underlying matrix is a sparse matrix:
You can then use that however you wish, with some of the tools in DifferentialEquations.jl.
You probably will face performance issues in your first implementation: I see from your current code that you create a lot of temporary objects, where the cost of computation stops being the actual numerics and starts being dominated from reserving and freeing memory (i.e. you have a lot of “allocation” overhead in your code). That is something we can address as well, probably making the code 10x or 100x faster, but we can discuss this after you have your first prototype working.