are tough! Usually the sum cannot be computed exactly.

Instead of your series we may take

a:=1/(i^4+j^4);

S:=Sum(Sum(a,i=1..infinity),j=1..infinity);


(this differs by yours by constant which is easy to compute).

value(S); does not work, as expected

Unfortunately Maple cannot approximate its sum dirrectly.

evalf(S); #does not work (unexpected!)

But the inner sum can be computed exactly, ==> 

evalf(  sum(sum(a,i=1..infinity),j=1..infinity) );

So, the series seems to be convergent, but this is not a proof.

( It is not difficult to see that se series is convergent using  a <= 1/i^2 * 1/j^2 ).

Another way for convergence:

S is convergent iff J < oo, where

J:=Int(Int(a,i=1..infinity),j=1..infinity):

Unfortunately J cannot be computed directly. 

value(J); # does not work

We must use:

int(a,i=1..infinity) assuming j::posint:
int(%,j=1..infinity);

so, J<oo and hence S<oo.

To compute a partial sum and verify numerically the sum:

eval(  'evalf'( 'add'('add'(a,i=1..N),j=1..N)),  N=1000);
     0.7617308252

(close to the value obtained earlier).

with(IntegrationTools):
f:=t^2*exp(-t)/ln(t):
J1:=Change(Int(f,t=0..1-a),t=(1-a)*x,x):
J2:=Change(Int(f,t=1+a..infinity), t=(1+a)/x,x):
simplify(limit(combine(J1+J2),a=0));

evalf(%);

      2.153547692

f := (i,j) -> `if`(j=1, i^2, i^3):
Matrix(10,2,f);

The system is too complicated to be solved. Use fsolve instead, it works.

Do the same for "Cello" (and conclude).

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.