The integral is computed by ScalarPotential  (if it exists, i.e. the form is exact).
Here is an example and a calling procedure (it works for more than 2 variables).

int1:=proc(F::list, X::list(name))
  uses VC=VectorCalculus;
  VC:-SetCoordinates('cartesian'[X[]]);
  VC:-ScalarPotential(VC:-VectorField(F))
end:

int1([2*x*y^7, 7*x^2*y^6], [x,y]);
      x^2*y^7

If you want inert integrals, it's easy to modify ScalarPotential.

It seems that the multiplication of float matrices with dimensions > 2000 crashes Maple 2021.

restart:
n:=2001:  # n:=2000 is ok
A:=Matrix(n, (i,j)->(i+2*j)/1000., datatype=float[8]):
A^2: #crash

Ortho.mw

A basis always exists, but the coefficients are probably huge and many primes are needed.

Symbolic computation with floats is inherently problematic.

f:=piecewise(x<0, exp(1), exp(1.)):
diff(f,x);

             piecewise(x = 0., Float(undefined), 0.)

restart:
a:=1:
before := {anames('user')};
#                         {a}
b:=2:
after := {anames('user')} minus eval(before,1) minus {'before'};
#                         {b}

You want

series(f/f1, a=0, 7);
expand(%);   # optional

but you have missing terms and other mistakes in your expected result. 

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