Solving an algebraic-differential equation

I cannot find the relevant part in the documentation. Maybe you have a quick solution to my problem.

Background

Sometimes you want to propagate in time some system of ODEs, but not all the variables are needed. Think of a matrix Schrödinger equation. Let us say, a sum of squares of unknown functions is desired. One may still solve ODEs of interest, then build the desired observable. But this is also an ideal case for an algebraic-differential solver. One just need to propagate everything, but having only the observable as an output.

My problem is the following messages

NDSolve::mconly: For the method IDA, only machine real code is available. Unable to continue with complex values or beyond floating-point exceptions. NDSolve::icfail: Unable to find initial conditions that satisfy the residual function within specified tolerances. Try giving initial conditions for both values and derivatives of the functions. 

Minimal working example

h[t_]:={{1,1,5},{1,2,I},{5,-I,3}} o={{0,0,5},{0,0,0},{5,0,0}} dim=3; vars=y[#]&/@Range[dim]; varsT=y[#][t]&/@Range[dim]; eqs=MapThread[Equal,{-I D[varsT,t],h[t].varsT}]; y0={0,1,0}; ics=MapThread[Equal,{y[#][0]&/@Range[dim],y0}; r=NDSolve[Join[eqs,ics],vars,{t,0,10}]//First; 

As the result we get some oscillating functions, an observable can be easily computed and plotted

ψ=vars/.r oT=Sum[Conjugate[ψ[[i]][t]]o[[i,j]]ψ[[j]][t],{i,dim},{j,dim}]; Plot[oT,{t,0,10}] 

I want to avoid post-processing to evaluate the observable oT.

Minimal not working example

Instead of evaluating the observable after solving ODEs, let us propagate it along with the unknowns. We add one additional algebraic equation and the respective initial condition

eqA={τ[t]==Sum[Conjugate[y[i][t]]o[[i,j]]y[j][t],{i,dim},{j,dim}]}; icA={τ[0]==0}; 

However, this fails with the error messages above

NDSolve[Join[eqs,eqA,ics,icA],τ,{t,0,10}] 

What can be done here? I am not even sure the initial condition icA for the auxiliary variable τ is needed.