NDSolve mixing many scalar and vector equations

I have a set of scalar equations in many unknowns, which I want to combine with a vector equation inside an NDSolve. The equations are a mix of differential and algebraic equations.

A set of scalar equations:

eqns = {-a[0][t] + a[1][t] == 0, a[1]'[t] == - fx[1][t] + fy[1][t] - a[1][t],
        a[2]'[t] == - fx[2][t] + fy[2][t] - a[2][t], a[2][t] + a[3][t] == 0,
        0 == x[0][t], a[1][t] == - x[0][t] + x[2][t], a[2][t] == - x[1][t] + x[3][t], 
        a[3][t] == x[1][t] + x[3][t], 1 == y[0][t], a[1][t] == - y[0][t] + y[2][t], 
        a[2][t] == - y[1][t] + y[3][t], a[3][t] == -y[2][t] + y[3][t]};

Along with a vector equation, defined using a function fCalc.

feqn = {{fx[0][t], fy[0][t]}, {fx[1][t], fy[1][t]}, {fx[2][t], fy[2][t]}, {fx[3][t], fy[3][t]}} == 
      fCalc[{{x[0][t], y[0][t]}, {x[1][t], y[1][t]}, {x[2][t], y[2][t]}, {x[3][t], y[3][t]}}];

My actual function fCalc is complicated, but it only evaluates for numerical input. It takes a list of {x,y} points and returns a list of {x,y} points, for instance:

fCalc[pts_ /; MatrixQ[pts, NumericQ]] := pts

Some initial conditions and the variables used:

initcs = {a[0][0] == 0, a[1][0] == 0, a[2][0] == 0, a[3][0] == 0};
vars=Flatten[Table[{a[j], x[j], y[j], fx[j], fy[j]}, {j, 0, 3}]];

Trying to solve this directly gives an error, as Mathematica tries to evaluate this vector function before starting, and refuses to start as it doesn't see enough equations for the unknowns.

NDSolve[{eqns, feqn, initcs},vars, {t, 0, 1}];
NDSolve::underdet: There are more dependent variables, than equations, so the system is underdetermined.

Removing the numerical requirement the system is well-defined:

fCalc2[pts_] := pts
AbsoluteTiming[NDSolve[{eqns, feqn /. fCalc -> fCalc2, initcs}, vars, {t, 0, 1}];]

So my question is, how to I trick Mathematica into waiting to evaluate the NDSolve until after giving values to the variables, and without having to repeatedly evaluate the fCalc function. There is an answer for solving for a vector system at this question, but I can't write my whole system in terms of a derivative of a vector as in that case.

My actual system is much more complicated (see the initial version of this question if you want to see), but I definitely need the numerical restriction on fCalc.