No problem, no worries.
I had an idea, which should work in principle. I'm pretty sure it fails because I ran out of memory. It actually almost completely freezes my laptop, and I have to force-quit Mathematica and several other applications to get back to where I can use it normally. That makes me think that the problem with your code is similar, except that your code fails semi-gracefully (no error reported, but does not freeze computer). Often running out of memory results in an error message, but I recall times when it did not.
I have only 16GB, but if you have more, the idea is to recast the system in terms of a vector variable cc comprising all your c[i, j] variables:
ca = CoefficientArrays[Sys[[;; Length@Bc, 2]], Through[Bc[t]]]
VectorQ[#["NonzeroValues"], NumericQ] & /@ ca (* check there are no variables in the coefficients *)
(* {True, True, True} *)
ivs = Sys[[Length@Bc + 1 ;;, 2]]; (* initial values *)
VectorQ[ivs, NumericQ] (* check there are no variables in the values *)
(* True *)
MemoryConstrained[
NDSolve[{
cc'[t] == Fold[#1.cc[t] + #2 &, Reverse@ca], (* reconstruct polynomial RHS of the ode system (Horner's method) *)
cc[0] == ivs},
cc, {t, tMin, tMax; 0.001}(*Method->"StiffnessSwitching", MaxSteps -> Infinity*)],
4*^9, (* put your memory limit here *)
"ran out of memory"]
(* "ran out of memory" *)
Should you be able to get a solution this way, the output value of the interpolation function will be a vector whose components correspond to the variables in Bc in the same order.
