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!

I had a brief look at your large system. You seem to want to plot a function of 2 variables x and t.
Here is what I believe is a simple example of the situation as far as I understand it:

restart;
sys:=diff(a(t),t,t)=a(t)+b(t)+a(t)*b(t),
     diff(b(t),t,t)=sin(t) - a(t);
res:=dsolve({sys,a(0)=0,D(a)(0)=0,b(0)=0,D(b)(0)=0},numeric);
plots:-odeplot(res,[t,a(t)],0..10); #No problem, but no x
plots:-odeplot(res,[t,b(t)],0..4.2);
#However you want to plot something like a(t)*x+b(t)*x^2
#If you are satisfied with an animation with x as the animation parameter then you can do:
p:=(x,T)->plots:-odeplot(res,[t,a(t)*x+b(t)*x^2],0..T);
plots:-animate(p,[x,7],x=0..1);
################################
#Alternatively you can use output=listprocedure, which will also enable you to use plot3d on the expression
#a(t)*x+b(t)*x^2
res2:=dsolve({sys,a(0)=0,D(a)(0)=0,b(0)=0,D(b)(0)=0},numeric,output=listprocedure);
A,B:=op(subs(res2,[a(t),b(t)])); #Numerical procedures  for a and b
plots:-odeplot(res2,[t,a(t)],0..10); #odeplot can still be used
plot(A,0..10); #This just gives  you the same as above
#Plotting a(t)*x+b(t)*x^2 for a fixed t = t0:
t0:=3.4567:
plot(A(t0)*x+B(t0)*x^2,x=0..1);
#Finally a 3d plot:
plot3d(A(t)*x+B(t)*x^2,x=0..1,t=0..7);
################################
#The singularity in your system may be quite real, i.e. it is actually there and not just due to numerical problems. It is to be expected because of the terms quadratic in the a[i]'s.
In the simple example above, try (after adding maxfun=0):
plots:-odeplot(res,[t,b(t)],0..62);
#A nonnumeric and even simpler example of a solution which has a finite interval of definition:
dsolve({diff(x(t),t)= 1+x(t)^2,x(0)=0});
#output x(t) = tan(t)

#################################
Trying to avoid the problems with your large system (described by PatrickT) I crudely set all floats to 3, keeping their signs, however.
SYS1:=evalindets(SYS,float,3*signum): #floats set to 3
M1:=dsolve(SYS1 union inits,numeric, output=listprocedure):
assign(seq(AA[i],i=1..25)=op(subs(M1,[seq(a[i](t),i=1..25)]))): #Cannot use A as that is used in your worksheet as a matrix.
plots:-odeplot(M1,[t,a[1](t)],0..2.1); #Works as in the simple example
#mm is 5
#3d plot:
plot3d(x*add(AA[m](t)*x^(m-1),m=1..mm),x=0..1,t=0..2,axes=boxed);