SolveTools:-SemiAlgebraic({c*x^7-6*x^3+6=0}, parameters=[c]);

piecewise(c < 0, [[x = RootOf(_Z^7*c-6*_Z^3+6, index = real[1])]], c = 0, [[x = 1]], c < RootOf(823543*_Z^3-1492992, 1.188 .. 1.25), [[x = RootOf(_Z^7*c-6*_Z^3+6, index = real[1])], [x = RootOf(_Z^7*c-6*_Z^3+6, index = real[2])], [x = RootOf(_Z^7*c-6*_Z^3+6, index = real[3])]], c = RootOf(823543*_Z^3-1492992, 1.188 .. 1.25), [[x = RootOf(50421*_Z^5*RootOf(823543*_Z^3-1492992, 1.188 .. 1.25)^2+62208*_Z^6+74088*_Z^4*RootOf(823543*_Z^3-1492992, 1.188 .. 1.25)-117649*_Z^2*RootOf(823543*_Z^3-1492992, 1.188 .. 1.25)^2+108864*_Z^3-172872*_Z*RootOf(823543*_Z^3-1492992, 1.188 .. 1.25)-254016, index = real[1])], [x = RootOf(50421*_Z^5*RootOf(823543*_Z^3-1492992, 1.188 .. 1.25)^2+62208*_Z^6+74088*_Z^4*RootOf(823543*_Z^3-1492992, 1.188 .. 1.25)-117649*_Z^2*RootOf(823543*_Z^3-1492992, 1.188 .. 1.25)^2+108864*_Z^3-172872*_Z*RootOf(823543*_Z^3-1492992, 1.188 .. 1.25)-254016, index = real[2])]], RootOf(823543*_Z^3-1492992, 1.188 .. 1.25) < c, [[x = RootOf(_Z^7*c-6*_Z^3+6, index = real[1])]])

P.S. It could be instructive to solve it by hand; it is not difficult, writing the equation as  c = f(x).

Edit. If you want only positive roots, use:

SolveTools:-SemiAlgebraic({c*x^7-6*x^3+6=0,x>0}, parameters=[c]);

Replace:

plot(x, x=0..1);
with
print(plot(x, x=0..1));

(this is because of how printlevel works, see ?printlevel).
 

seq(
evalf(Int(abs(J), [T[2] = T2start .. T2end, x[5] = 0 .. a], epsilon=1e-4)),
a=[10^3,10^4,10^5]
);

So, your integral seems to be divergent

To obtain directly the wanted result just use:

f:=10000/(1+30762*0.478^x)+5:
maximize(diff(f,x), location);

    1845.361367, {[{x = 14.00001597}, 1845.361367]}

Now let me explain why those complex numbers appeared.
The graph of the derivative is

Due to inherent round-offs (because you did not use rationals and solve instead of fsolve) your assigned M could be a bit bigger than the actual maximum of f'. So, the equation f'(x)=M will have only complex roots. To get rid of the imaginary part you may take Re(a), where a is the complex root.

Your matrix T is a standard test for eigenvector computation. It is symmetric but has “repeated” eigenvalues.

Maple can find the exact eigenvalues/eigenvectors in this case.

V,Q:=Eigenvectors(T);

But the test is designed for approximate computations (floating point).
So, use:
V1,Q1 := Eigenvectors(evalf(T));

(Then take Re to get rid of the small imaginary parts due to round-offs).

Note that in both cases, the matrix Q (or Q1)  containing the eigenvectors is not orthogonal (even if such a matrix exists for symmetric matrices). Eigenvectors() finds it only for float symmetric matrices with distinct eigenvalues.
If  orthogonality is needed then use GramSchmidt.
 

Download ode-ex.mw

Denoting F(n) = f(0)+...+f(n) :

rsolve({f(n+1)=F(n), f(n)=F(n)-F(n-1), F(0)=1}, {f(n),F(n)});

(f(n) for n>0 of course)

p:=ln((1+x)*y) + exp((x^2)*(y^2)) - x - cos(x):
n:=40: Order:=n+1: interface(rtablesize=infinity):
Y:=convert(solve(series(p,x),y),polynom):
Vector[column](n, k ->  (D@@k)(y)(0)=eval(diff(Y,x$k),x=0));

(D(y))(0) = 0
((D@@2)(y))(0) = -2
((D@@3)(y))(0) = -2
((D@@4)(y))(0) = 55
((D@@5)(y))(0) = 96
((D@@6)(y))(0) = -3951
((D@@7)(y))(0) = -14986
((D@@8)(y))(0) = 599292
((D@@9)(y))(0) = 3756040
((D@@10)(y))(0) = -150583195
((D@@11)(y))(0) = -1451365992
((D@@12)(y))(0) = 56622924821
((D@@13)(y))(0) = 787829738580
((D@@14)(y))(0) = -29495018843850
((D@@15)(y))(0) = -570264628995262
((D@@16)(y))(0) = 20174284393035491
((D@@17)(y))(0) = 528406101043111104
((D@@18)(y))(0) = -17368560725263003371
((D@@19)(y))(0) = -607502621004862751222
((D@@20)(y))(0) = 18165803871525900350124
((D@@21)(y))(0) = 844845810439069562531252
((D@@22)(y))(0) = -22331822067408858287791607
((D@@23)(y))(0) = -1391102518495130466232499136
((D@@24)(y))(0) = 31137730503244680996639906697
((D@@25)(y))(0) = 2661297665528307975243093795336
((D@@26)(y))(0) = -46909537305923595629791853669346
((D@@27)(y))(0) = -5812207204576583472301320482485946
((D@@28)(y))(0) = 69478963661245363181701637264539327
((D@@29)(y))(0) = 14236970345056860768024681386350762560
((D@@30)(y))(0) = -71299313756801875315745344773429106215
((D@@31)(y))(0) = -38343398748183919497538363485168620149282
((D@@32)(y))(0) = -148767292026687176665730760219959487334644
((D@@33)(y))(0) = 110608691261039295372096825657213070145581600
((D@@34)(y))(0) = 1784186834273881762868359210688981525024869565
((D@@35)(y))(0) = -327307575749373732034049337639072936118525057752
((D@@36)(y))(0) = -11596412918978194063904912726445275723035187765539
((D@@37)(y))(0) = 899618765531255422817716130916175428672344561662044
((D@@38)(y))(0) = 67838929477109472349986391474357426166573238064107894
((D@@39)(y))(0) = -1466744386949476812393980478836845534302116051221857142
((D@@40)(y))(0) = -389148909149217184540519589123668475583410103807659331061

I prefer:

deq := FunctionAdvisor(DE, hypergeom([a, b], [c], z))[2, 1]:
Deq:=parse(convert(deq,string)):
dsolve(Deq, f(z));

T:=u-> degree( eval(u,diff=1), [f(x),g(x)]):
select( u -> T(u)<2, p);

In evalf[20], a = .1  is computed with Digits=10.
You should use subs(a=1/10) or Digits:=20

# ex := 2.*ca*g+ka*ma+ka*meqc+keqc*ma-1.*ma*meqc*(ca*keqc+2.*g*ka)/(ca*ma+ca*meqc+g*ma): #TableauRouth[3,1]
ex := 
(2*ca^3*g*keqc*ma+2*ca^3*g*keqc*meqc+4*ca^2*g^2*ka*ma+4*ca^2*g^2*ka*meqc+2*ca^2*g^2*keqc*ma+ca^2*keqc^2*ma^2+4*ca*g^3*ka*ma+2*ca*g*ka^2*ma^2+4*ca*g*ka^2*ma*meqc+2*ca*g*ka^2*meqc^2+ca*g*ka*keqc*ma^2-3*ca*g*ka*keqc*ma*meqc+ca*g*keqc^2*ma^2+2*g^2*ka^2*ma^2-2*g^2*ka^2*ma*meqc+g^2*ka*keqc*ma^2)*(2*ca^2*g*ma+2*ca^2*g*meqc+2*ca*g^2*ma+ca*ka*ma^2+2*ca*ka*ma*meqc+ca*ka*meqc^2+ca*keqc*ma^2+g*ka*ma^2-g*ka*ma*meqc+g*keqc*ma^2): #[4,1]
Explore( g=solve(ex>0,g), parameters=[ca=0..5.0, ka=0..5.0, keqc=0..5.0, ma=0..5.0, meqc=0..5.0], echoexpression=false );

heunB is 1 for r=0 and continuous. The integral is divergent due to 1/r.

The inequalities are polynomial (after eliminating the denominators: a/b>0 <==> a*b>0).
So they can be solved with SemiAlgebraic or even solve, but a HUGE and useless result should be expected.