You can just use  simplify  command. For a more substantial simplification seems to be  specified numerical values ​​to all constants.

simplify(dsolve({ODE,bcs}));

SeqInColumn:=proc(S)

local i;

for i to nargs-1 do

print(args[i], ``);

od;

args[nargs];

end:

Your example:

B := b=2:

C := c=3:

A := B, C:

SeqInColumn(A);

                b=2,
                 c=3

Your example:

restart;

eq:=f(x)=1+x*f(x)^2;

f(x):=a+b*x+c*x^2+d*x^3;

expand(eq);

sys:={seq(coeff(lhs(eq), x, k)=coeff(rhs(eq), x, k), k=0..3)};

solve(sys);

Or

member(5, {seq(13^k mod 100, k=1..1000)});

                                   false

The benefit a set compared to a list that automatically repetitive elements removed from a set and it is sorted in ascending order.

Try

 Hits := identify(d);

for i from 1 to 501 do if type(Hits[i], 'float') = false then  print(Hits[i]) end if end do;

The law of sines:

AB/sin(C)=BC/sin(A)=AC/sin(B)=2*R=6

We have  sin(C)=AB/6=sqrt(2)/2,  sin(A)=BC/6=sqrt(2)/3

Two cases are possible:  C=Pi/4  or  C=3*Pi/4

Uniquely  A=arcsin(sqrt(2)/3)  because  A<C

Obviously that  cos(B)=cos(Pi - (A+C))=-cos(A+C)  

Next we use Maple and the law of cosines:

c1:=expand(-cos(Pi/4+arcsin(sqrt(2)/3)));

c2:=expand(-cos(3*Pi/4+arcsin(sqrt(2)/3)));

AB:=3*sqrt(2): BC:=2*sqrt(2):

radnormal(sqrt(AB^2+BC^2-2*AB*BC*c1));  # first answer

radnormal(sqrt(AB^2+BC^2-2*AB*BC*c2));  # second answer

Addition - visualization of both solutions:

expand(2*sin(x+(1/4)*Pi));
``*coeff(%, sin(x))*sin(x)+``*coeff(%, cos(x))*cos(x);

Addition:  In Standard M 16 the same result can be obtained  simpler

 expand(``*2*sin(x+(1/4)*Pi));

subs(2=`2`, expand(2*(sin(x+Pi/4))));

Likely this can only be done numerically:

f := BesselJ(1, t)*(z*BesselJ(0, t)*BesselJ(1, z)-t*BesselJ(0, z)*BesselJ(1, t))/((1+10*t)*(z^2-t^2)):

F:=unapply(Int(f, t=0..infinity), z);

evalf([F(0), F(1), F(2)]);

The original formulation of the problem needs to be clarified. The  maximize  command  finds the maximum value of the function in the speified region. If the region is not specified, then in  the whole domain. But the function  F := (x,y)-> (-x*y+12)*x*y/(2*x+2*y)  in  R^2  can take any real value from  -infinity  to  infinity . This is easily seen if we consider the restriction of this function on the line  y=x :  F(x,x)=x*(12-x^2)/4 . Properties of the plot of a cubic polynomial are well know:

plot(x*(12-x^2)/4, x=-5..5);

Task to become more interesting, we look for the maximum value of the function  F(x,y)  in the region  x>=0, y>=0 . Using Maple it's easy to find the critical points of  F (x, y)  in this region:

F:=(x,y)->x*y*(12-x*y)/2/(x+y):

solve({D[1](F)(x,y)=0, D[2](F)(x,y)=0, x>=0, y>=0}); 

                                      {x = 2, y = 2}

But the proof that at this point the function  F(x,y)  has the maximum value in the region  x>=0, y>=0  there is not very simple. We carry it in 2 stages:

1. Extend the function at the origin  F(0,0)=0 , and show that the function is continuous in the region  x>=0, y>=0. Enough to prove the continuity in the point  (0,0) . This  follows from the simple estimates. Denote  d=sqrt(x^2+y^2) . We have

x*y<=(x^2+y^2)/2=d^2/2 ,   x+y>=sqrt(x^2+y^2)=d ,   x*y*(12-x*y)/2/(x+y)<=12*d^2/(2*d)=6*d

2. Next we estimate the value of the function  F(x,y)  at an arbitrary line  x+y=C , C>0 :

x*y*(12-x*y)/2/(x+y)<=36/(2*C)=18/C  , because  maximize(t*(12-t))=36

Thus, the maximum value of a continuous function  F(x,y)  is achieved in the triangle  0<=y<=6-x ,  0<=x<=6  and it is now clear that it will  F(2,2)=4

The problem is obvious - find the sum in closed form:

sum((2*k-1)^2, k=1..n);

factor(%);

We introduce the notation:  V1=<1,0,0>,  V2=<-1,1,0>,  V3=<0,0,1>  ,  V=<3,0,1>. The vectors  V1 , V2 , V3  are the basis vectors in R^3. Any linear mapping is uniquely determined by its values ​​on the basis vectors. 

So you will need to do two steps:

1. Expand the vector  V  by these basis vectors, ie find numbers  x1, x2, x3  that  V=x1*V1+x2*V2+x3*V3 

2. Use the property of the linear mapping T, ie  T(V)=x1*T(V1)+x2*T(V2)+x3*T(V3)

restart:

N :=unapply( N0-(1/2)*sqrt(2)*sqrt(Pi*Kc/d)*Sum(erfc((1/2)*(L*n+x)/sqrt(d*t)) + erfc((1/2)*((n+1)*L-x)/sqrt(d*t)), n = 0 .. infinity), x);

N0:=0.2: L:=0.25: Kc:=2*10^(-12): t:=360000: d:=2.010619298*10^(-10):

plot(N,0..0.25,axes=box, labels=[x, N]);

if  V  is a eigenvector for the Matrix  A  then  A.V=lambda*V ,  lambda  is corresponding eigenvalue. So

A:=<<3,-2,-1,2,0>|<11,4,-8,2,7>|<0,0,2,0,0>|<3,3,-4,3,3>|<-8,4,5,-4,-1>>:

V:=<2,0,-1,0,1>:  V1:=A.V:

print(` A`=A,`   V`=V,`   V1`=V1);

``;

lambda:=V1[1]/V[1];

It is easy to prove the special case of your identity, if we take  n0=n,  n1=0  (n>=1) :

sum((-1)^(i-1)*factorial(n)/(i^2*factorial(i-1)*factorial(n-i)), i = 1 .. n);

sum(1/i, i = 1 .. n);

                                               Psi(n+1)+gamma

                                               Psi(n+1)+gamma

Now apply this identity first for n = n0, then for n = n1, and then subtract the second equality of the first one. 

Finally replace the summation index in the resulting equality.