Your system is polynomial. I have substituted g = G*K.
The system reduces to find the roots of a huge polynomial in I1 of degree 9, with 8 parameters.
There will be 9 solutions. A few hundred pages will be enough for the result (if you arrange them nicely).

restart;
F:=[b*R + b*S - d*S - bcd*I1*S - bcd*I2*S - (e*g*S)/K - (g*I1*S)/K - ( g*I2*S)/K - (g*R*S)/K - (g*S^2)/K,
-d*e - (e^2*g)/K - (e*g*I1)/K - (e*g*I2)/K - (e*g*R)/K + bcd*I1*S +  bcd*I2*S - (e*g*S)/K - e*scd1,
-d*I1 - (e*g*I1)/K - (g*I1^2)/K - (g*I1*I2)/K - (g*I1*R)/K - ( g*I1*S)/K + e*scd1 - e*f*scd1 - I1*scd2,
-acd*I2 - d*I2 - (e*g*I2)/K - (g*I1*I2)/K - (g*I2^2)/K - ( g*I2*R)/K - (g*I2*S)/K + e*f*scd1,
-d*R - (e*g*R)/K - (g*I1*R)/K - (g*I2*R)/K - (g*R^2)/K - (g*R*S)/K +  I1*scd2]:
P:=eval(F, g = G*K):
indets(P);
#      {G, I1, I2, R, S, acd, b, bcd, d, e, f, scd1, scd2}
V:=[S,e,I1,I2,R];
#                     V := [S, e, I1, I2, R]
eliminate(P,{S,e,R}):
sys:=[%[2][]]:
map(degree, sys, [I1,I2]);
#                             [5, 5]
eliminate(sys,{I2}):
polyI1:=%[2][]:
degree(polyI1,I1);
#                               9
nops(polyI1);
#                               10
length(polyI1);
#                            507655
convert(polyI1,string): length(%);
#                             432373
indets(polyI1);
#             {G, I1, acd, b, bcd, d, f, scd1, scd2}

Maple is not able to recognize a particular 0 expression. This is not uncommon in any CAS. In our case:

ex := -sqrt((-8*x - 16)*exp(x/2) + x^2 + 4*x + 16*exp(x) + 4) + 4*exp(x/2) - x - 2:

is(simplify(ex)=0) assuming x::real;
        FAIL
simplify(eval(ex, x=2*t)) assuming t::real;  #workaround
        0

I don't think this is a true bug. The designer simply did not anticipate that a user would choose numerical labels. Probably such numerical labels should not be accepted.

You cannot prescribe all three points of tangency, because this implies 6 conditions and the ellipse has only 5 parameters.
But it is possible (in general) to fix two tangency points and the third will be computed.

You must check the definitions when comparing the results from different CAS-es.
Zeta(0, 3, 0.5);
      8.414398322

restart
ode:=diff(y(x),x) = 2*(2*y(x)-x)/(x+y(x)): ic:=y(0)=2:
mu:=DEtools:-intfactor(ode):
sol:=dsolve([mu*ode,ic], [exact], implicit):
simplify(map(exp,sol));

To avoid complex solutions, let's assume x>0,y>0.

1. ode_original can be factored. This will explain the difficulties:

(2*x^(5/2)-3*y(x)^(5/3))  *  (2*(diff(y(x), x))*x-3*y(x))  /  (6*x^(5/2)*y(x)^(5/3)) = 0

The first factor does not depend on y'   and ==>  y = a * x^(3/2) for a = ...  [the singular solution]
The second factor ==> your linear ode ==> y = C x^(3/2), so, it contains the previous singular solution.

2. Now, about the dsolve solution without simplification. All the listed solutions except the last are of the form y = b * x^(3/2) for some b.
The last solution can be written as  z^(-3/2) + z = - C1,  where z = y / x^(3/2).
It rezults that z = C2,  so,  y = C2 * x^(3/2). Hence, it contains indeed the general solution.

In packages such as VectorCalculus or DifferentialGeometry  the functions are usually supposed to be "sufficiently smooth".
For a "pathological" function it's better to use the definition.
The directional derivative of F in the direction of the unit vector V is:
Limit((F(X+t*V)-F(X))/t, t=0);
In our case:

F:=(x,y) -> min(abs(x),abs(y)):
(F(0-t/sqrt(2), 0+t/sqrt(2)) - F(0,0))/t;
limit(%, t=0);

        undefined

It is a bug. If mysol is given as   lhs-rhs=0  or lhs-rhs, it works.

odetest((lhs-rhs)(mysol) = 0, [ode,ic]);
      [0,0]
 

simpSum:=proc(S::specfunc(Sum))
local o1:=op(1,S), eq:=op(2,S), n,a,b,c;
if not type(eq, name=range) then return S fi;
n:=lhs(eq); a,b:=op(rhs(eq));
c:=max(select(type, [solve(o1,n)], integer));
if c>=a and c <= b then
  sum(o1, n=a..c-1) + Sum(o1, n=c+1..b)
else S fi;
end:

simpSum(Sum(n*x^n, n = 0 .. 20));

S:=Sum(n*(n-4)*x^(n^2), n = 0 .. infinity);
simpSum(S);

S := Sum(epsilon^k*v(k,t)/k!, k = 0 .. infinity);

S1 := diff(S, epsilon);

                       
 S1ok:=simpSum(S1);                   

limit(%, epsilon=0);

The result is correct. diff(x^n,x)  is 0 for n=0, not for n=1.
diff(Sum(x^n, n = 0 .. infinity), x) is correct too. Note that diff(Sum(f,...), x)  simply applies diff to f.
Of course, it is possible to simplify by eliminating the 0 term, but unfortunately it has to be done by hand (if you want).

S:=diff(Sum(x^n, n = 0 .. infinity),x);
eval(S, 'n = 0 .. infinity' = 'n = 1 .. infinity');

(In this case you may simply use eval(S, 0=1).)

I think that setting interface(typesetting=extended) solves the problem, but I cannot check this in Maple 2015.
But another problem is that the continued fraction appears as 1/( 1/... + 1). To change this, use inert operators g := x -> 1/(1 %+ x)

1 + something   is not  1.  It corresponds to Not(something).

with(GraphTheory):
G := RandomGraphs:-RandomGraph(8, .5):
u:=3;
DrawGraph(G);
add(map2(Degree, G, Neighbors(G,u)));