Definitely the error is due to your "..." (not shown).

You may have used it incorrectly, but Groebner solves at once your problem.

sys:=O1*O3*f7-f1, O1*O3*f1-f7, O1*O3*f8-f2, O1*O3*f2-f8, O1*O3*f5-f3, O1*O3*f3-f5, O1*O3*f6-f4, O1*O3*f4-f6, O1*O4*f8-f1, O1*O4*f1-f8, O1*O4*f7-f2, O1*O4*f2-f7, O1*O4*f6-f3, O1*O4*f3-f6, O1*O4*f5-f4, O1*O4*f4-f5, O1*O2*O3*f6-f1, O1*O2*O3*f1-f6, O1*O2*O3*f5-f2, O1*O2*O3*f2-f5, O1*O2*O3*f8-f3, O1*O2*O3*f3-f8, O1*O2*O3*f7-f4, O1*O2*O3*f4-f7, O1*O2*O4*f5-f1, O1*O2*O4*f1-f5, O1*O2*O4*f6-f2, O1*O2*O4*f2-f6, O1*O2*O4*f7-f3, O1*O2*O4*f3-f7, O1*O2*O4*f8-f4, O1*O2*O4*f4-f8:

oo:=O1^2-1,O2^2-1,O3^2-1,O4^2-1:

Groebner[Basis]([sys,oo], plex(f1, f2, f3, f4, f5, f6, f7, f8,  O1, O2, O3, O4)):
select(has,%,[seq(f||i,i=1..8)]);

You may solve eqs wrt all the variables:

eqs := {x - a*y,y - a*x}:
solve(eqs);

n:=3:
Matrix(combinat[permute]( [-1$n, 1$n], n) );

seq(x||i,i=1..200);

So, you have a homogeneous linear system of ODEs with constant coefficients.
If you know a fundamental system for the solutions:

exp(p*t), exp(q*t), ...

where p,q, ... are distinct complex numbers
then the characteristic polynomial in variable y is

(y - p)(y - q) ...

[if there are multiple roots, the situation is a bit more complicated]

Note that c(t) = 0 cannot appear in a fundamental system
(because of the linear independence).

indets(eq,anyfunc(identical(t)));

If you must use fsolve, try to restrict the domain.

fsolve({f, g}, {a, c}, a = 10000 .. infinity, c = 0 .. infinity);

Edit.
1. You may use
plots[implicitplot]([f,g], a=-20..20,c=-20..20,color=[red,blue]);
to localize the roots
2. The system seems to be intentionally constructed to have "non-intuitive" solutions.

Yes, obviously for the ics case Maple "forgot" the BesselJ term .

M := diff(T(r), r, r)+(diff(T(r), r))/r+u*T(r)+P*(r^4+r^2) = 0; # u includes the constant

M0 := {M, D(T)(0)=0}:
s:=dsolve(M);

s0:=dsolve(M0);


s0general:=T(r)=BesselJ(0, sqrt(u)*r)*_C2-P*(r^4*u^2+r^2*u^2-16*r^2*u-4*u+64)/u^3;

odetest(s, M);
                               0
odetest(s0, M0);
                              {0}
odetest(s0general, M0);
                              {0}

Edit. If another ics is added e.g.
M01 := {M, D(T)(0)=0, T(0)=a};
then BesselJ appears!

B:=Groebner[Basis]([x-v, y-v^2], plex(u,v,x,y)):
remove(has,B,[u,v]);

Here is an ad-hoc procedure.

flatt:=proc(L::{list,set})
  local a,n; a:=L; n:=-1;
  while nops(a)<>n do
    n:=nops(a);
    a:=map(x->`if`(type(x,{list,set}),op(x),x),a);
  od;
  op(a)
end:

flatt([1,2,3, {4,5,6}]);

ex:=-Omega*a*sqrt(2)*sqrt(-Omega^2*a^2-2*k*m+sqrt(Omega^2*a^2*(Omega^2*a^2+4*k*m)))/(-Omega^2*a^2+sqrt(Omega^2*a^2*(Omega^2*a^2+4*k*m))):

radnormal(ex^2);
      -1

So, ex = I or -I
Note that both values may appear, depending on the parameters: e.g. changing a to -a  ==> ex to -ex
ex = I for  e.g. Omega=1, a=2, k=-1, m=1.

Edit. As Markiyan has pointed out, OP says Omega, a, k, m are >0. (I did not see this because the horizontal scroll was absent in my browser). But in this case it is easy to see that ex = -I. In fact, the denominator is obviously >0,  so

must be <0 (because ex^2=-1).

It follows that the argument of ex is - Pi/2 ==> ex = - I.

The solution for your ODEs is very oscillatory.
Maple can easily find a series solution, but it will be almost useless.
I think that the solution itself would be useless if found.

Please look at the following similar simple example:

Digits:=30: Order:=10:
a:=10^6:
s1:=dsolve({diff(y(t),t)=cos(a*t)*y(t), y(0)=1}, y(t), series):
s2:=dsolve({diff(y(t),t)=cos(a*t)*y(t), y(0)=1}, y(t)):
ss1:=convert(rhs(s1),polynom);

ss2:=rhs(s2); #exact

plot(ss2, t=0..1);  #exact

plot( ss1, t=0..1); #series