Maple will solve inequalities with parameters in them, but only if the expressions are rational functions.  So, the sqrt in your problem makes things tricker.

Here is an example of that where your expression has been multiplied by its conjugate:

expr:=0 < 2*b^2*(10*a*K+3*K*b-sqrt((K+(2*a*K+1)/b)^2-4*K/b)*b+5)*(10*a*K+3*K*b+sqrt((K+(2*a*K+1)/b)^2-4*K/b)*b+5);
solve( expand~(expr), {a}, parametric);

@tomleslie Playing around, Digits:=30; might be enough but Digits:=25 is not enough.

Just a note:  while solve may sometimes give a bare RootOf() as an answer, the usage pattern of allvalues( RootOf(equation, variable) );  is almost never better than calling solve(equations, variable); unless you know exactly what you are doing.

I think the command you are looking for is msolve.

Are you just trying to test ideal membership or are you trying to do something more esoteric than that?

You could start at the help page for repetition statemetns in Maple https://www.maplesoft.com/support/help/maple/view.aspx?path=do if you need more help than that, you may need to include more details about what you are trying to do.

Joe keeps a GitHub project of the Emacs Maple Mode.  Last updated two years ago I think. https://github.com/JoeRiel/maplev

Check your constraints as a first step. It looks like Maple's Optimization routines can't find any values which satisfy those (even with integer and non-negative assumptions removed).

@Zeineb Even if you install the third-party DirectSearch package, I doubt it will be able to solve your n=100 problem (10,000 equations!).  I suspect you need to reformulate to something that can be handled by a specialty numerical library.

@mehdibaghaee I would have to see your full problem to say for sure, but if you are doing assignments in a loop, you could instead do evaluations in the loop -- but in that case you should use map[inplace](eval, M, {a=0}) (e.g.) to avoid copying at matrix at each evaluation.  Another option might to use rtable_eval(Mm__11, inplace) to force evaluation of the assigned symbols to modify the Matrix.

@tomleslie If you read more carefully, you'll see it's only four unknowns.  The @bogus5331 has just not provided us with those values explicitly (you can see that only the four unknowns in the fsolve command appear in the pasted ascii 2D-output)

PolynomialTools:-Splits will get you part of the way there, but then you're faced with the problem of identifying real linear factors, and combining the complex ones.   I don't, right off, know an out of the box solution that gives a purely real factorization.

Something to be aware of.  For each value of epsilon, your system has (generically) has four solutions.  fsolve by default only returns one of the solutions to a system. See with:

sol := solve({eq3, eq4, eq5, eq6}, {d[0], d[1], d[2], d[3]});
sols := allvalues(sol):
evalf(eval([sols], epsilon = -1/10));

It is possible to call Matlab solvers, but in fact Matlab mostly uses all the same numerical numerical solvers that Maple does ( ATLAS / MLK / LAPACK ).   The key to getting at the speed in Maple is making sure your matrix is stored as a hardware floating point (datatype=float[8]) array so that you avoid expensive copy and conversion overhead when you call linear algebra routines.

@AmusingYeti creating a function f as a procedure and then plotting it as a function call  plot(f(x), x=0..5) is a classic anti-pattern in Maple, it works OK in many cases, but you are almost always better off plotting the procedure directly:

plot(f, x=0..5);
plot(2*f, x=0..5);

etc