You can set the environmental variable UseHardwareFloats to false.
Simple example:

restart;
Digits:=8:
UseHardwareFloats:=false:
LinearAlgebra:-GaussianElimination(Matrix([[1.2345,6.789],[11.1213,-9876]]));

See
?LinearAlgebra,details
at the bottom there is a link to a worksheet with examples: Examples/LA_NAG

The method of Frobenius involves using infinite series, so you may try
ode:= x*diff(y(x),x,x)+2*diff(y(x), x)+16*x*y(x) = 0;
dsolve(ode,y(x),formal_solution);
value(%);

Since this is homework I shall try not to give you all the details.
Since you are asked to find the equilibrium solutions (when H < A it should say) it is convenient to start out by defining the right hand side as a function f.

restart;
f:=y->A*y-B*y^2-H*y;
ode:=diff(y(t),t) = f(y(t));
res1:=dsolve({ode,y(0)=y0});
#Solve f(y)=0 w.r.t. y to find the equilibria.
param:={A=1,B=1,H=0.2,y0=2};
#You could now plot eval(rhs(res1),param) just to see how the population develops with constant harvesting.
#To introduce intermittent harvesting you could use the following function of t:
h:=add((-1)^i*Heaviside(t-3*i),i=0..3);
plot(h,t=0..10,thickness=3,discont=true);
#Now replace H by h*H:
ode2:=eval(ode,H=h*H);
#Solve and plot as before
#You can convert the result (res2) involving Heaviside into piecewise by:
convert(res2,piecewise) assuming t>0;

Does this do it?

plot3d([2*t^3-t^4,(t+2)*cos(theta)-2,(t+2)*sin(theta)],t=0..2,theta=0..2*Pi,scaling=constrained);
plots:-spacecurve([t,-2,0],t=0..1.7,color=red,thickness=3);
plots:-display(%%,%);

There is a multiplication sign missing. If you insert that you get

(0.1250000000e-3*(.6000000000*Unit('m')+.6510936027*(Unit('m'))*(.2929921212*Unit('m'))/(Unit('m')*((1/2)*Pi-26.38779995))))*Pi^2*Unit('m')^2*((1/2)*Pi-26.38779995);

and

simplify(%);
results in -0.1813470337e-1*Unit('m'^3).

eq:=diff(y(x),x,x)+y(x)^2=x^4+2;
sol:=dsolve({eq,y(0)=0,y(1)=1},numeric);
plots:-odeplot(sol,[x,y(x)],0..1);
#Second equation
eq:=diff(y(x),x,x)+y(x)^2=-b*y(x)^2*diff(y(x),x);
sol:=dsolve({eq,y(0)=1,D(y)(0)=0},numeric,parameters=[b]);
sol(parameters=[.1]);
plots:-odeplot(sol,[x,y(x)],0..1);

You could use fsolve:

plot(tan(sqrt(lambda)) - sqrt(lambda),lambda=0..200);
fsolve(tan(sqrt(lambda)) = sqrt(lambda), lambda = .1 .. 30);
fsolve(tan(sqrt(lambda)) = sqrt(lambda), lambda = 50..70);

inequal({x^2+y^2<= 5^2, 2^2 < x^2+y^2, arctan(y,x) < -3*Pi*(1/4) or 3*Pi*(1/4) < arctan(y,x)}, x=-5..5,y=-5..5,axiscoordinates = polar);

Modifying my answer to your earlier question:

restart;
with(LinearAlgebra):
K:=Matrix(3,symbol=b);
V:=Vector(3,symbol=a);
Vt:=apply~(V,t);
eq1:=a[1](t)*a[2](t)^2+a[2](t)*a[3](t)*a[1]+a[3](t)^3;
eq2:=a[1](t)^3*a[2](t)^2+a[2](t)*a[3](t)+a[2](t)^3;
eq3:=a[1](t)*a[3](t)^2+a[3](t)^2*a[1](t)^3+a[2](t);
W:=<eq1,eq2,eq3>;
K.Vt=W; #Need to solve for K, i.e. the b's:
Equate(K.Vt,W);
res:=solve(%,{seq(seq(b[i,j],j=1..3),i=1..3)}); #Using solve
K1:=subs(res,K);
simplify(K1.Vt);
#Any K1 does the job, but since you want the matrix K1 to be invertible you cannot just set the free variables (the remaining b's) to zero. So you could try:
K2:=eval(K1,{seq(seq(b[i,j]=i+2*j,j=1..3),i=1..3)});
simplify(K2.Vt);
Determinant(K2);

Will the following work for your system?
In this example there are 2 dependent variables x and y. I use RootFinding:-NextZero on their product.

restart;
sys:=diff(x(t),t)=sin(y(t))+sin(t), diff(y(t),t)=-sin(x(t))-y(t);
res:=dsolve({sys,x(0)=0,y(0)=Pi/2},numeric,output=listprocedure);
X,Y:=op(subs(res,[x(t),y(t)]));
plots:-odeplot(res,[[t,x(t)],[t,y(t)]],0..20);
Txy:=NULL:
t1:=0:
to 5 do t1:=RootFinding:-NextZero(t->X(t)*Y(t),t1); Txy:=Txy,t1 end do;
Vt:=Vector([Txy]);
Vx:=X~(Vt);
Vy:=Y~(Vt);
M:=Matrix([Vt,Vx,Vy]);

eq:=y(n+2) + 7/10 *y(n+1) + 1/10* y(n) = 1;
#Using rsolve
rsolve({eq,y(0) = 0,y (1) = 1},y(n));
#Using ztrans:
ztrans(eq,n,z);
solve(%,{ztrans(y(n),n,z)});
eval(%,{y(0) = 0,y (1) = 1});
invztrans(op(%),z,n);

Here is a way:

A:=Int(sin(x)^3/x^3,x=0..infinity);
combine(A);
IntegrationTools:-Parts(%,-sin(3*x)+3*sin(x));
IntegrationTools:-Parts(%,-3*cos(3*x)+3*cos(x));
#And now write the result as a sum of two integrals. On one make a very simple change of variable.
#Compare with
value(A);

restart;
dC1 := diff(C1(t), t) = -a1*C1(t)+b1*C2(t);
dC2 := diff(C2(t), t) = a1*C1(t)-(b1+a2)*C2(t)+b2*C3(t);
dC3 := diff(C3(t), t) = a2*C2(t)-(b2+a3)*C3(t)+b3*O1(t);
dO1 := diff(O1(t), t) = a3*C3(t)-b3*O1(t);
syst := dC1, dC2, dC3, dO1;
ics := C1(0) = 1, C2(0) = 0, C3(0) = 0, C1(t)+C2(t)+C3(t)+O1(t) = 1;
#Your odes and the initial conditions imply that C1(t)+C2(t)+C3(t)+O1(t) = 1, so I would leave out O1:
simplify(dC1+dC2+dC3+dO1);
#Using method=laplace you get a rather messy looking answer which suggests using numerical solution:
sol3:=dsolve({dC1, dC2, eval(dC3,O1(t)=1-(C1(t)+C2(t)+C3(t))), ics[1..3]},{C1(t), C2(t), C3(t)},method=laplace);
#For a numerical solution you need concrete values for your rate constants:
params:={a1=.1,a2=.2,a3=.3,b1=.123,b2=.345,b3=.456}; #Arbitrarily chosen
solN:=dsolve(eval({dC1, dC2, eval(dC3,O1(t)=1-(C1(t)+C2(t)+C3(t))), ics[1..3]},params),numeric);
plots:-odeplot(solN,[[t,C1(t)],[t,C2(t)],[t,C3(t)],[t,1-C1(t)-C2(t)-C3(t)]],0..50,legend=[C1,C2,C3,O1]);
#However, the laplace solution is just fine:
L:=eval(sol3,params);
plot(subs(L,[C1(t),C2(t),C3(t),1-C1(t)-C2(t)-C3(t)]),t=0..50,legend=[C1,C2,C3,O1]);

You need to be very cautious when moving among square roots. Better not to simplify than do an unwarranted one:
u:=sqrt(x*y)/(x*y);
simplify(u); #No change
w:=sqrt(z)/z;
#Outputs w:=1/sqrt(z) so in fact you get the result you want by doing
eval(w,z=x*y);
#That caution is needed is well known, but here is an example:
u1:=sqrt(x*y)/(x*sqrt(y));
simplify(u1); #No change
u2:=simplify(u1,symbolic);
eval(u1,{x=-1,y=-1});
eval(u2,{x=-1,y=-1});
#So is it reasonable to expect Maple (or any other program) to see that in u there appears the common product x*y in numerator and denominator? Or rather, is it worth the effort to program for occurrences like that?

For your first question:

L:=[[1,3], [5,8], [9,2], [5,12]];
#Either
map2(op,2,L);
#or
op~(2,L);

For the second question you could use the sequence version suggested by Konstantin, or the more elaborate
L2:=[4,2,-1,9,5,6,2,0,14,7];
L2e:=ListTools:-Enumerate(L2);
map(x->`if`(x[1]::even,x[2],NULL),L2e);

That one is not particularly elegant, but more entertaining!