ode := diff(y(x), x) = sin(y(x))+1;
ic := y(0) = Pi;
Order := 40;
sol := dsolve([ode, ic], y(x), series);
ser1 := gfun[seriestorec](subs(sol, y(x)), a(n));
diffeq1 := gfun[rectodiffeq](ser1[1], a(n), y(x));
dsol1 := dsolve(diffeq1, y(x));

This should help:

indets(cS(a)(b)(c), specfunc(anything, cS));

Maple does not resolve types in the 0th operand.

Replace
Int(unapply(eq7,t),-25..25,method=_d01ajc)
with
Int(eq7,t=-25..25,method=_d01ajc)

Instead of ifactor use ifactors. Much easier to work with.

y = -1/2 .. 3, z = -3 .. 3, x = 0 .. 1;
u = 2*x*y, v = 2*x*z;
u = (minimize .. maximize)(2*x*y, x = 0 .. 1, y = -1/2 .. 3);
v = (minimize .. maximize)(2*x*z, x = 0 .. 1, z = -3 .. 3);
u = -1 .. 6, v = -6 .. 6, x = 0 .. 1;
uvEQ := collect(combine(subs(solve({u = 2*x*y, v =  2*x*z}, {y, z}), expand([eq1, eq2]*x^2*300))), [exp, sin, cos]);

Digits := 16;
PLT1:=plots[intersectplot](
op(map(E->map(evalf, surface(E, u = -1 .. 6, v = -6 .. 6, x = 0 .. 1)), uvEQ)),
maxlev = 5, maxtet=100000, grid = [31, 61, 11], thickness = 3, transparency = 0.3):
print(PLT1);

subs(solve({u = 2*x*y, v = 2*x*z}, {y, z}), [y, z, x]);
F1 := unapply(%, [u, v, x]);
plots[display](plottools[transform](F1)(PLT1), view=[ -1/2 .. 3,  -3 .. 3,  0 .. 1]);

If you'd like a plot of the solution you can use:

Digits := 16;
plots[intersectplot](
map(evalf, surface(eq1, y = -1/2 .. 3, z = -3 .. 3, x = 0 .. 1)),
map(evalf, surface(eq2, y = -1/2 .. 3, z = -3 .. 3, x = 0 .. 1)),
maxlev = 5, grid = [31, 61, 11], thickness = 3, transparency = 0.3);

You could use this, if you mean exact match.

ArrayTools:-IsZero(A) and LinearAlgebra:-Dimensions(A)=(4,4);

Also see ?ArrayTools:-IsEqual

expand(rel(n+1));

This works:

restart;
expr := -(r0+Delta_r)^2*(46*r0-41*Delta_r)*r0^5;
subsop(2=b,3=c,4=d, 1=a, expr);

Reasoning is obvious. subsop(1=a, expr) replaces the constant of multiplication. Operand 1 (always the constant) stops being a constant, so a new expression with 2 elements is formed: a*(rest of expression). There is no operand 3 and 4. Your subsop worked as intended.

Try this code, if you have the newer Maple version (looks like you're trying to use old Maple).

y:=sqrt((x^2*(x^2-4))/(x^2-1));
solve(evalc(Im(y)), x);
F := [Minimize = minimize, Maximize = maximize];
G := 'simplify(evalc(Re(y)))';
F(`assuming`([G], [-1 <= x, x <= 1]), x = -1 .. 1, location);
F(`assuming`([G], [2 <= x]), x = 2 .. infinity, location);
F(`assuming`([G], [x <= -2]), x = -infinity .. -2, location);

Simply submitting the relevant query to msolve will do.

{x = 26730899*a[2]+1162213*a[5]+507520*a[9]+808519*a[4]+757969*a[8]+252655*a[3]+12167*a[1]+529*a[7], y = 8046000599*a[2]+748465172*a[5]+3822524979*a[9]+8821208032*a[4]+588083951*a[8]+5934310109*a[3]+8346050986*a[1]+7538374844*a[7]+4385029649*a[6]};

{a[5] = 0, a[2] = 0, a[3] = 0, a[4] = 0, a[8] = 0, a[1] = 0, a[6] = 6880584, a[9] = 4174, a[7] = -4004513};

Do you know your expression is annihilated by P:

P := 1+256*u*(256*u^3-135)*lambda^4+(16384*u^3+1728)*lambda^3+1536*lambda^2*u^2+64*lambda*u;
`@`(evala, Normal, subs)(lambda = expr, P);
 

I would guess you are meant to supply irreducible forms. You can always try this:

(degree, factor)(evala(Norm(lambda-bSol[1])));
Digits:=150;
lcm(op(map(p->PolynomialTools:-MinimalPolynomial(evalf(p), x, 36), bSol)));

Or start with this:

PA := evala(Algfield(bSol)):
PA[3];
 

If you map objective via exp function the following works naturally:

Optimization:-Maximize(x[1]*x[2]*x[3], {0 <= x[1], 0 <= x[2], 0 <= x[3], x[1]+x[2] <= 1, x[1]+x[3] <= 2});

I did not anticipate this limitation of Maple's solve. Given these responses, ?solve, parametric is the better choice.
Example:

solve({0 < 4*t1+4*r1, t1^2-4*t2 = r1^2, 0 <= r1}, {t2, r1}, parametric, parameters = {t1});
solve({0 < 4*t1+4*r1, t1^2-4*t2 = r1^2, 0 <= r1}, {t2, t1}, parametric, parameters = {r1});