Example:

Max:=-infinity:

S:=[2, -5, 4, 0, 4, 1]:

for s in S do

if  s>Max  then  Max:=s  fi;

od:

Max;

                                                        4

To find the positive roots of your equations, I first put them to the reduced form of quadratic equation, and then find the positive roots of these equations using well-known formula, taking the plus sign in front of the square root:

eq1 := 1.6*10^(-7)*R*sin(t)-4.4*10^(-14)*R^2*cos(t)*sin(t)-1.6*10^(-14)*R^2*cos(t)^2+4.2*10^(-14)*R^2-1.3+2.1*10^(-9)*R*cos(t) = 0:

eq2 := 8.3*10^(-8)*R*sin(t)-1.2*10^(-13)*R^2*cos(t)*sin(t)-2.9*10^(-44)*R^2*cos(t)^2+7.1*10^(-14)*R^2-1.3+8.3*10^(-8)*R*cos(t) = 0:

eq3 := 8.3*10^(-8)*R*sin(t)-1.2*10^(-13)*R^2*cos(t)*sin(t)-2.2*10^(-44)*R^2*cos(t)^2+7.1*10^(-14)*R^2-1.3+8.3*10^(-8)*R*cos(t) = 0:

eq4 := 2.1*10^(-9)*R*sin(t)-4.4*10^(-14)*R^2*cos(t)*sin(t)+1.6*10^(-14)*R^2*cos(t)^2+2.6*10^(-14)*R^2-1.3+1.6*10^(-7)*R*cos(t) = 0:

L := [eq1, eq2, eq3, eq4]:  # List of original equations

L1 := map(collect, map(t->t/coeff(lhs(t), R^2), L), R):   # List of reduced equations

Sol := map(t->-(1/2)*coeff(lhs(t), R)+sqrt((1/4)*coeff(lhs(t), R)^2-coeff(lhs(t), R, 0)), L1):   # List of positive roots of this equations  

A := plot(Sol, t = 0 .. 2*Pi, color = [blue, red, red, green], thickness = 3, coords = polar):  #  Required curves in polar coordinates

B := plot(min(op(Sol)), t = 0 .. 2*Pi, color = yellow, coords = polar, filled = true): # The intersection of the regions

plots[display](A, B);  # All together

 Addition:   min(f(t), g(t))  -  the minimum function of  f(t)  and  g(t)

Obviously the equality  holds

x=<0.3603, 1.3279, 1.6882> = 0.3603*<1, 0, 1> + 1.3279*<0, 1, 1> 

From this equality it follows that  x  is in the rowspace  of the vectors  <1,0,1>  and  <0,1,1>

For the second and third graphics positive roots for any  t  almost are the same, so in the plot we see  not 4 but 3 graphics (blue graph covers the yellow graph):

eq1:=1.6*10^(-7)*R*sin(t)-4.4*10^(-14)*R^2*cos(t)*sin(t)-1.6*10^(-14)*R^2*cos(t)^2+4.2*10^(-14)*R^2-1.3+2.1*10^(-9)*R*cos(t)=0:

eq2:=8.3*10^(-8)*R*sin(t)-1.2*10^(-13)*R^2*cos(t)*sin(t)-2.9*10^(-44)*R^2*cos(t)^2+7.1*10^(-14)*R^2-1.3+8.3*10^(-8)*R*cos(t)=0:

eq3:=8.3*10^(-8)*R*sin(t)-1.2*10^(-13)*R^2*cos(t)*sin(t)-2.2*10^(-44)*R^2*cos(t)^2+7.1*10^(-14)*R^2-1.3+8.3*10^(-8)*R*cos(t)=0:

eq4:=2.1*10^(-9)*R*sin(t)-4.4*10^(-14)*R^2*cos(t)*sin(t)+1.6*10^(-14)*R^2*cos(t)^2+2.6*10^(-14)*R^2-1.3+1.6*10^(-7)*R*cos(t)=0:

L:=[eq1,eq2,eq3,eq4]:

L1:=map(collect,map(t->t/coeff(lhs(t),R^2), L),R):

plot(map(t->-coeff(lhs(t),R)/2+sqrt(coeff(lhs(t),R)^2/4-coeff(lhs(t),R,0)),L1), t=0..2*Pi, color=[red,blue,yellow,green], thickness=2);

 You can easily finish the remaining yourself.

restart;

S:=b*h^2/6: M:=10^7: fb:=M/S: fall:=2400: obj:=abs(fb-fall);

eval(obj, {b=10, h=50});  # Minimum

limit(eval(obj, {b=1/n,h=2/n}), n=infinity);  # Maximum

This works:

a:=t^2;

f:=unapply(a, t);

                                     a := t^2

                                  f := t -> t^2

Example:

f(3);

                 9

Definition of the function and its plotting for  alpha=2:

f := (x, y, alpha)->piecewise((-(1/2)*Pi <= x and x < 0 or 0 < x and x <= (1/2)*Pi) and -abs(x) <= y and y <= abs(x), (x^2+y^2)^alpha*arcsin(y/x),  x = 0 and y = 0, 0,  undefined);

plot3d(f(x, y, 2), x = -3 .. 3, y = -3 .. 3, numpoints = 25000, style = surface, axes = normal, orientation = [30, 70]);

Changing  the number  alpha  you will receive a variety of functions and their graphs.

restart;

with(plots):

B := 1.5:

u2 := r*cos(th): v2 := r*sin(th):

w1 := u1+I*v1: w2 := u2+I*v2:

z1 := evalc(w1^2): z2 := evalc(w2^2):

x1 := evalc(Re(z1)): x2 := evalc(Re(z2)):

y1 := evalc(Im(z1)): y2 := evalc(Im(z2)):

f1 := proc (theta) options operator, arrow; plot3d([-x1, -y1, v1], u1 = -6 .. 1, v1 = -B .. B, grid = [50, 50], orientation = [theta, 80], color = u1) end proc:

f2 := proc (theta) options operator, arrow; plot3d([-x2, -y2, v2], r = 1 .. 1, th = -Pi .. Pi, grid = [50, 50], orientation = [theta, 80], color = black, thickness=2) end proc:

display(seq(display(f1(10*k), f2(10*k)), k = -17 .. 18), insequence = true, axes = box, view = [-1 .. 1, -1 .. 1, -B .. B]);

eq:=diff(x(t),t,t) - 3*diff(x(t),t) + 2*x(t) = 2*t -3;

inc:=x(0)=2, D(x)(0)=4;

dsolve({eq, inc});

In radians Maple calculates immediately:

cos(Pi/6),  tan(Pi/4);

                                        1/2*3^(1/2),  1

If the angles are given in degrees, it is necessary to convert to radians:

cos(convert(30*degrees, radians)), tan(convert(45*degrees, radians));

                                        1/2*3^(1/2),  1

A, B, C := <x , 2 , -1>, <0, y , 1>, <3 , 4 , -2>:

AB, BC, CA := B-A, C-B, A-C:

d := u->LinearAlgebra[Norm](u,2,conjugate=false):

p := (u,v)->LinearAlgebra[DotProduct](u,v,conjugate=false):

Sol := [evalf(RealDomain[solve]({p(CA,AB)=0,d(CA)=d(AB)})),

evalf(RealDomain[solve]({p(AB,BC)=0,d(AB)=d(BC)})),

evalf(RealDomain[solve]({p(BC,CA)=0,d(BC)=d(CA)}))];

V:=map2(eval, [A,B,C], Sol):

T:=seq(plottools[polygon](map(convert,V[i],list)), i=1..4):

plots[display](T, axes=normal, scaling=constrained, orientation=[50,70]);

Edit: at first I did not notice that only needs to be AB = CA. Therefore, from the solutions should be left only the first two.

Example if  2 circles in the same plane:

eq1:=<0.8*cos(t),0.8*sin(t),0>:

eq2:=<2*cos(t),2*sin(t),0>:

C1:=plots[spacecurve](eq1, t=0..2*Pi, color=blue, thickness=3):

C2:=plots[spacecurve](eq2, t=0..2*Pi, color=blue, thickness=3):

S:=plot3d(s*eq1+(1-s)*eq2, t=0..2*Pi, s=0..1, style=surface, color=yellow):

plots[display](map(plottools[rotate],map(plottools[translate],[C1,C2,S], 3,4,1),Pi/3,0,Pi/6), axes=normal, orientation=[-15,60]);

 This method is suitable not only for the concentric circles and even if they lie in different planes:

eq1:=<0.8*cos(t)+0.5,0.8*sin(t),1>:

eq2:=<2*cos(t),2*sin(t),0>:

C1:=plots[spacecurve](eq1, t=0..2*Pi, color=blue, thickness=3):

C2:=plots[spacecurve](eq2, t=0..2*Pi, color=blue, thickness=3):

S:=plot3d(s*eq1+(1-s)*eq2, t=0..2*Pi, s=0..1, style=surface, color=yellow):

plots[display](map(plottools[rotate],map(plottools[translate],[C1,C2,S], 3,4,1),Pi/3,0,Pi/6), axes=normal, orientation=[120,45], scaling=constrained);

I think that to solve your problem it is convenient to use  Optimization  package:

restart;

de1 := diff(V(t), t) = V(t)-(1/3)*V(t)^3-W(t)+Ie:

de2 := diff(W(t), t) = 0.8e-1*(V(t)+.7-.8*W(t)):

sys:={de1, de2}:

sol:=dsolve({op(eval(sys,Ie=0.5)), V(0)=0, W(0)=1}, numeric, output=listprocedure);

Optimization[NLPSolve](rhs(sol[2]), 0..200, maximize, method=branchandbound );

Optimization[NLPSolve](rhs(sol[2]), 0..200);

beta := 1: lambda := 5: nu := 1: Delta := 1: 

sys := {C(0) = 0, In(0) = .2, diff(C(u), u) = beta*In(u)*s(u)-C(u)/nu, diff(In(u), u) = C(u)/nu-In(u)/lambda+Delta*(diff(diff(In(u), u), u)), diff(s(u), u) = -beta*In(u)*s(u), s(0) = .8, (D(In))(0) = .2}: 

dsys := dsolve(sys, numeric, [In(u), s(u), C(u)]): 

plots[odeplot](dsys, [[u,In(u)], [u,s(u)], [u,C(u)]], u = 0 .. 0.5, thickness=2, axes = boxed, legend=[In(u),s(u),C(u)]);

Because part of the plot is not included in the domain of the function  x, we extend the plot by  line segments using piecewise  command.

restart;

Eq1 := x = phi*(theta[m]-1/x)/theta[m]+(1-phi)*(1/2):

Sol := solve(Eq1, x):

x := unapply(Sol[1], phi):

p1 := plot(piecewise(`and`(theta[m] > .9, theta[m] < 1), 5*theta[m]-4.5, theta[m] > 1, x(0)), theta[m] = .9 .. 14, color = red, legend = "a=0 (extreme bimodal)"):

a := eval(x(.1), theta[m] = 1.55):

p2 := plot(piecewise(`and`(theta[m] > .9, theta[m] < 1.55), a*(theta[m]-.9)/(.65), theta[m] > 1.55, x(.1)), theta[m] = 1 .. 14, color = green, linestyle = 3, legend = "a=0.1"):

b := eval(x(.5), theta[m] = 3.56):

p3 := plot(piecewise(`and`(theta[m] > 3.55, theta[m] < 3.56), b*(theta[m]-3.55)/(0.1e-1), theta[m] > 3.56, x(.5)), theta[m] = 3.55 .. 14, color = blue, linestyle = 6, legend = "a=0.5"):

p4 := plot(piecewise(`and`(theta[m] > 3.9, theta[m] < 4), 5*theta[m]-20, theta[m] > 3.9, x(1)), theta[m] = 3.9 .. 14, color = black, linestyle = 2, legend = "a=1 (uniform)"):

P1 := plot([map(proc (t) options operator, arrow; [t, eval(x(0), theta[m] = t)] end proc, [1, 2, 8, 12])], style = point, color = red, symbol = cross, symbolsize = 17):

P2 := plot([map(proc (t) options operator, arrow; [t, eval(x(.1), theta[m] = t)] end proc, [1.55, 2.5, 3.5, 8, 12])], style = point, color = green, symbol = diagonalcross, symbolsize = 17):

P3 := plot([map(proc (t) options operator, arrow; [t, eval(x(.5), theta[m] = t)] end proc, [3.56, 5.5, 6, 8, 12])], style = point, color = blue, symbol = asterisk, symbolsize = 17):

P4 := plot([map(proc (t) options operator, arrow; [t, eval(x(1), theta[m] = t)] end proc, [4, 5, 8, 12])], style = point, color = black, symbol = box, symbolsize = 15):

plots[display]({P1, P2, P3, P4, p1, p2, p3, p4}, axes = boxed, view = [0 .. 14, 0 .. 1]);