Since your expression is a square trinomial with respect to x, it is sufficient to solve the system  {f(a,b)>0, discrim(Delta,x)<0} 

Example:

restart;
Delta:=f(a,b)*x^2+g(a,b)*x+h(a,b):
Sys:={f(a,b)>0, discrim(Delta,x)<0};
f(a,b):=1-a:  g(a,b):=a+b:  h(a,b):=a-b:
solve(Sys);

See help on  ?procname(recursion)

We need to help Maple a bit:

f:=(5*x^2+3*(x+exp(x))^(1/3)+exp(x)*(2*x^2+3*x))/x/(x+exp(x))^(1/3);
IntegrationTools:-Expand(int(f, x));
simplify(%);  # The result
simplify(diff(%, x));  # Checking

Addition. In fact, the second integral  int(1/(sqrt(x+exp(x))), x)  is more hard. Mathematica fails with it also.

For  fsolve  command, the number of equations must equal the number of unknowns. I took  V=1 :

restart;
omega1 := 1.562:
omega2 := 2.449:
omega3 := 3.325:
y1 := c1*sin(omega1*t+phi1)+c2*sin(omega2*t+phi2)+c3*sin(omega3*t+phi3): 
y2 := .1019*c1*sin(omega1*t+phi1)+.75*c2*sin(omega2*t+phi2)+.4608*c3*sin(omega3*t+phi3):
y3 := .407*c1*sin(omega1*t+phi1)+(0*c2)*sin(omega2*t+phi2)+1.844*c3*sin(omega3*t+phi3): 
eq1 := subs(t = 0, y1) = 0:
eq2 := subs(t = 0, y2) = 0: 
eq3 := subs(t = 0, y3) = 0: 
eq4 := subs(t = 0, diff(y1, t)) = V:
eq5 := subs(t = 0, diff(y2, t)) = 0:
eq6 := subs(t = 0, diff(y3, t)) = 0:
V:=1: 
eqs := {eq1, eq2, eq3, eq4, eq5, eq6}: 
vars := {c1, c2, c3, phi1, phi2, phi3}: 
fsolve(eqs, vars);

    {c1 = -0.8218023979, c2 = 0.0001356814559, c3 = 0.08520994334, phi1 = -9.424777961, phi2 =    3.141592654, phi3 = -3.141592654}

You can save V variable without assigning any value to it. In this case, you must use  solve  command and you will get several solutions:

restart;
omega1 := 1.562:
omega2 := 2.449:
omega3 := 3.325:
y1 := c1*sin(omega1*t+phi1)+c2*sin(omega2*t+phi2)+c3*sin(omega3*t+phi3): 
y2 := .1019*c1*sin(omega1*t+phi1)+.75*c2*sin(omega2*t+phi2)+.4608*c3*sin(omega3*t+phi3):
y3 := .407*c1*sin(omega1*t+phi1)+(0*c2)*sin(omega2*t+phi2)+1.844*c3*sin(omega3*t+phi3): 
eq1 := subs(t = 0, y1) = 0:
eq2 := subs(t = 0, y2) = 0: 
eq3 := subs(t = 0, y3) = 0: 
eq4 := subs(t = 0, diff(y1, t)) = V:
eq5 := subs(t = 0, diff(y2, t)) = 0:
eq6 := subs(t = 0, diff(y3, t)) = 0: 
eqs := {eq1, eq2, eq3, eq4, eq5, eq6}: 
vars := {c1, c2, c3, phi1, phi2, phi3}: 
solve(eqs, vars);

restart;
evalc((2-I*X)^4);
Re(%) assuming X::real;

or

restart;
expand(Re((2-I*X)^4)) assuming X::real;

or

restart;
evalc(Re((2-I*X)^4));

Probably the last method is the shortest  (evalc  assumes that all parameters are real numbers).

sort(D^4 + D^3 + D^2 +D + C, D,  descending);

For a solution in  Optimization  package it is useful to use the  initialpoint  option, which can be obtained from the plot:

restart;
plot3d(`if`(y >= -(1/13)*x-85/13 and y <= -(25/14)*x+26/7 and y <= (23/12)*x+89/12, 4*x*y-2*y^2+4*x+y+4, NULL), x=-10..10,y=-10..10, numpoints=40000, axes=normal);
Optimization[Maximize](4*x*y-2*y^2+4*x+y+4,{y >= -(1/13)*x-85/13, y <= -(25/14)*x+26/7, y <= (23/12)*x+89/12}, initialpoint=[x=-5,y=-5]);
Optimization[Minimize](4*x*y-2*y^2+4*x+y+4,{y >= -(1/13)*x-85/13, y <= -(25/14)*x+26/7, y <= (23/12)*x+89/12}, initialpoint=[x=3,y=-5]);

P:=x^8 + 2*x^4 + 1; 
Q:=x^5 + 4*x^3 + x^2 + 3*x +3;
3*P+4*Q mod 5;
P-Q mod 5;
expand(P*Q) mod 5;
R:=rem(P, Q, x, 'q') mod 5;  # remainder
q;  # quotient
is(P=expand(q*Q+R) mod 5);  # Check

Edit.
 

Use a semicolon or a colon after each assignment:

equations:={b[1](2)*a[1](2, 1) = 1/2, b[1](2)*a[1](2, 1)^2 = 1/3, b[1](2)*a[2](2, 1) = 1/2, b[1](2)*a[2](2, 1)^2 = 1/3, b[2](2)*a[1](2, 1) = 1/2, b[2](2)*a[2](2, 1) = 1/2, b[1](2)*(a[3](2, 1)+a[3](2, 2))*a[1](2, 1) = 1/3, b[1](2)*(a[3](2, 1)+a[3](2, 2))*a[2](2, 1) = 1/3, b[1](2)*a[2](2, 1)*a[1](2, 1) = 1/3, b[1](1)*(a[3](1, 1)+a[3](1, 2))+b[1](2)*(a[3](2, 1)+a[3](2, 2)) = 1/2, b[1](1)*(a[3](1, 1)+a[3](1, 2))^2+b[1](2)*(a[3](2, 1)+a[3](2, 2))^2 = 1/3, b[2](1)*(a[3](1, 1)+a[3](1, 2))+b[2](2)*(a[3](2, 1)+a[3](2, 2)) = 1/2, b[1](1)+b[1](2) = 1, b[2](1)+b[2](2) = 1};

variables:={a[1](2, 1), a[2](2, 1), a[3](1, 1), a[3](1, 2), a[3](2, 1), a[3](2, 2), b[1](1), b[1](2), b[2](1), b[2](2), b[3](1), b[3](2)};

solve(equations,variables);

If you want to build several curves defined by explicit equations in the same plot, use square brackets (a list of curves). Here is an example of plotting all 4 basic trigonometric functions:

restart;
plot([sin(x), cos(x), tan(x), cot(x)], x=-Pi/2..5*Pi/2, y=-4..4, color=[red, blue, green, gold], size=[800,500], scaling=constrained,  discont);
 

You can easily do this one line without repeating:

restart;
f(x):=exp(-x+3.):
T__4:=unapply(convert(taylor(f(x), x=7, 5), polynom), x);  # T__4 as an operator
T__4(3);

or (without unapply):

restart;
f(x):=exp(-x+3.):
T__4:=convert(taylor(f(x), x=7, 5), polynom);  # T__4 as an expression
eval(T__4, x=3);

To calculate the values of the derivative at individual points use  unapply  command or  D  operator:

KK := unapply(diff(x^2, x), x);
KK(2);  # or
D(x-> x^2)(2);

We just write the variable  s  as an index:

sys :=
  diff(x[s](t),t) = y[s](t),
  diff(y[s](t),t) + x[s](t) = 0;

ic := x[s](0) = a(s),
      y[s](0) = b(s);
dsolve({sys, ic});

restart:    
with(IntegrationTools):
M:=Int(lambda(s)*diff(u(s),s,s,s), s);
for n from 1 to 3 do
Parts(M, `if`(n=1, lambda(s), diff(lambda(s), s$(n-1))));  M:=%;
od:
M;

If you specify parameter values, then it is easy to find all the solutions numerically.

Example:

restart;
eq1 := (1/1152)*(81*a1*P^4*b^4*C1^3+162*a1*P^4*b^4*C1*C2^2+648*a1*P^4*b^4*C1*C3^2+432*a1*P^4*b^4*C1*C4^2+864*a1*P^4*b^4*C2*C3*C4+18*a2*P^4*a^2*b^2*C1^3+90*a2*P^4*a^2*b^2*C1*C2^2+90*a2*P^4*a^2*b^2*C1*C3^2+288*a2*P^4*a^2*b^2*C1*C4^2+81*a3*P^4*a^4*C1^3+648*a3*P^4*a^4*C1*C2^2+162*a3*P^4*a^4*C1*C3^2+432*a3*P^4*a^4*C1*C4^2+864*a3*P^4*a^4*C2*C3*C4+288*d1*P^4*b^4*C1+576*d2*P^4*a^2*b^2*C1+288*d3*P^4*a^4*C1+1152*d4*P^4*a^2*b^2*C1+4096*Tc*a^3*b^3*C4)/(b^3*a^3) = 0:

eq2 := (1/1152)*(162*a1*P^4*b^4*C1^2*C2+864*a1*P^4*b^4*C1*C3*C4+81*a1*P^4*b^4*C2^3+432*a1*P^4*b^4*C2*C3^2+648*a1*P^4*b^4*C2*C4^2+90*a2*P^4*a^2*b^2*C1^2*C2+72*a2*P^4*a^2*b^2*C2^3+612*a2*P^4*a^2*b^2*C2*C3^2+360*a2*P^4*a^2*b^2*C2*C4^2+648*a3*P^4*a^4*C1^2*C2+864*a3*P^4*a^4*C1*C3*C4+1296*a3*P^4*a^4*C2^3+432*a3*P^4*a^4*C2*C3^2+2592*a3*P^4*a^4*C2*C4^2+288*d1*P^4*b^4*C2+2304*d2*P^4*a^2*b^2*C2+4608*d3*P^4*a^4*C2+4608*d4*P^4*a^2*b^2*C2-4096*Tc*a^3*b^3*C3)/(b^3*a^3) = 0:

eq3 := (1/1152)*(648*a1*P^4*b^4*C1^2*C3+864*a1*P^4*b^4*C1*C2*C4+432*a1*P^4*b^4*C2^2*C3+1296*a1*P^4*b^4*C3^3+2592*a1*P^4*b^4*C3*C4^2+90*a2*P^4*a^2*b^2*C1^2*C3+612*a2*P^4*a^2*b^2*C2^2*C3+72*a2*P^4*a^2*b^2*C3^3+360*a2*P^4*a^2*b^2*C3*C4^2+162*a3*P^4*a^4*C1^2*C3+864*a3*P^4*a^4*C1*C2*C4+432*a3*P^4*a^4*C2^2*C3+81*a3*P^4*a^4*C3^3+648*a3*P^4*a^4*C3*C4^2+4608*d1*P^4*b^4*C3+2304*d2*P^4*a^2*b^2*C3+288*d3*P^4*a^4*C3+4608*d4*P^4*a^2*b^2*C3-4096*Tc*a^3*b^3*C2)/(b^3*a^3) = 0:

eq4 := (1/144)*(54*a1*P^4*b^4*C1^2*C4+108*a1*P^4*b^4*C1*C2*C3+81*a1*P^4*b^4*C2^2*C4+324*a1*P^4*b^4*C3^2*C4+162*a1*P^4*b^4*C4^3+36*a2*P^4*a^2*b^2*C1^2*C4+45*a2*P^4*a^2*b^2*C2^2*C4+45*a2*P^4*a^2*b^2*C3^2*C4+36*a2*P^4*a^2*b^2*C4^3+54*a3*P^4*a^4*C1^2*C4+108*a3*P^4*a^4*C1*C2*C3+324*a3*P^4*a^4*C2^2*C4+81*a3*P^4*a^4*C3^2*C4+162*a3*P^4*a^4*C4^3+576*d1*P^4*b^4*C4+1152*d2*P^4*a^2*b^2*C4+576*d3*P^4*a^4*C4+2304*d4*P^4*a^2*b^2*C4+512*Tc*a^3*b^3*C1)/(b^3*a^3) = 0:

S:=indets([eq1,eq2,eq3,eq4]) minus {C1,C2,C3,C4};  # The list of the parameters

assign(S=~{$1..nops(S)});  # Assigning numeric values to parameters

S;  # Values of the parameters

evalf([solve({eq1, eq2, eq3, eq4}, {C1, C2, C3, C4})]);  # All the numerical solutions