If we look closely at the marks on the t-axis, we see that all the graphs are constructed in a logarithmic scale. I did not understand why this happened, because I did not find the appropriate settings in your document. By default Maple uses a uniform scale on both axes. I just restarted your document and got other graphics that do not cause questions:

Download schechter_guo_v2_new.mw

In your code there are two "if", but only one "end if" .  Also what means  epsilon^() ?

The last initial condition  limit(diff((diff(f(x), x))/x, x), x = 0) = 0  follows from two conditions  D(f)(0)=0  and  (D@@2)(f)(0)=0 . But then the equation is obtained overdetermined (5 conditions instead of 4)  and Maple shows an error. The equation can be solved with the condition  D(f)(0)=0  and not very high accuracy:

restart;
ode := x^3*(diff(f(x), x, x, x, x))+alpha*(x^4*(diff(f(x), x, x))+x^3*(diff(f(x), x, x))-x^2*(diff(f(x), x)))-2*x^2*(diff(f(x), x, x, x))+3*x*(diff(f(x), x, x))-3*(diff(f(x), x))+R*x*f(x)^2-R*x^2*(diff(f(x), x, x))*(diff(f(x), x))-3*R*x*f(x)*(diff(f(x), x, x))+3*R*f(x)*(diff(f(x), x))+x^2*R*f(x)*(diff(f(x), x, x, x))-M^2*(x^3*(diff(f(x), x, x, x))-x^2*(diff(f(x), x))) = 0:
ics := f(0) = 0, f(1) = 1, D(f)(1) = 0, D(f)(0)=0:
sol:=dsolve(eval({ode, ics},[alpha=1,R=2,M=3]), method=bvp[midrich],numeric, maxmesh=1000, abserr=10^(-3));
plots:-odeplot(sol, [x,f(x)], x=0..1, color=red); 

Edit.

f:=(R,h,l)->[(1/3)*Pi*R^2*h, Pi*R*l, Pi*R^2+Pi*R*l]:

Examples of use:

f(3,4,5);
map(t->f(t[]), L);
                          

You do not need any loops. See this toy example:

A:=Matrix([[["BrakePower", "BrakePower"], [2.356, 2.356], [2.749, 2.749], [3.142, 3.142], [3.534, 3.534], [3.927, 3.927], [4.32, 4.32]], [["BrakePower", "S2"], [2.356, .303], [2.749, .271], [3.142, .256], [3.534, .249], [3.927, .244], [4.32, .241]], [["BrakePower", "S4"], [2.356, .256], [2.749, .225], [3.142, .211], [3.534, .205], [3.927, .2], [4.32, .197]]]):
y := a*x^2+b*x+c:
CurveFitting[LeastSquares](convert(A[1,2..-1], list), x, curve = y);
 

For your "polynomial" you can do so:

norm(subs(1.25=125, aa), infinity);

Here is a more general method that is suitable for a "polynomial" with any number of terms with fractional exponents:

max(abs~(map(p->`if`(type(evalf(op(1,p)),numeric),op(1,p),1),[op(aa)]))[]);
 

Edit.

Should be  s:=s+i  instead of  s:=s+a[i]

Doubt_on_list_new.mw

Explanation:  i  is a member of the list  a , and  a[i]  is the member of the list  a standing at  i_th place.

f:=x->1/(1+x^2):
f1:=[[0,f(0)], [1/3,f(1/3)]]:
f2:=[[1/3,f(1/3)], [2/3,f(2/3)]]:
f3:=[[2/3,f(2/3)], [1,f(1)]]:
plot([f(x), f1, f2, f3], x=0..1, color=[red,blue,green,black], view=[0..1,0..1], legend=['f','f1','f2','f3']);

In your question, you wrote  f(x)=1/1+x^2. But this is just  f(x)=1+x^2 . Therefore, I assumed that the parentheses are omitted.

Addition. If you want all functions to be plotted on the same interval  0 .. 1 , then do so:

f:=x->1/(1+x^2):
g:=(x1,y1,x2,y2)->y1+(y2-y1)*(x-x1)/(x2-x1):
f1:=g(0,f(0),1/3,f(1/3)):
f2:=g(1/3,f(1/3),2/3,f(2/3)):
f3:=g(2/3,f(2/3),1,f(1)):
plot([f(x), f1, f2, f3], x=0..1, color=[red,blue,green,black], view=[0..1,0..1.2], legend=['f','f1','f2','f3']);

Edit.
 

Should be

plot(eval([(exp((1/2)*x))*((1/2)*x+2*t*(1/3)), (exp((1/2)*x))*(1+t((-256)+240-108+27)/2^7+t(168-24+27)^2/2^7+(1/3)*t^3*((-16)+9)/2^6+(1/4)*t^4/2^7), (exp((1/2)*x))*(1+t(256*(-1.01)+240*(-1.01)^2+108*(-1.01)^3+27*(-1.01)^4)/2^7+t(168*(-1.01)^2+24*(-1.01)^3+27*(-1.01)^4)^2/2^7+(1/3)*t^3*(16*(-1.01)^3+9*1.01^4)/2^6+(1/4)*t^4*(-1.01)^4/2^7)], t=5), x = -4 .. 4);

Edit.

This can be done in several ways. The shortest one:

G := {A=exp1, B=exp2}:
assign(G);

A, B;

                                      exp1, exp2

If  alpha  occurs not in one, but in several places at once, then it is more convenient to use freeze  and  thaw commands.

Example:

restart;
expr:=-sin(alpha)*(sin(theta1)*cos(theta)-cos(theta1)*sin(theta))+2*cos(alpha)*sin(theta)*cos(theta):
combine(subs([sin(alpha)=freeze(sin(alpha)), cos(alpha)=freeze(cos(alpha))], expr)):
thaw(%);

                              sin(alpha)*sin(-theta1+theta)+cos(alpha)*sin(2*theta)

Obviously this is a bug. The problem (assuming all variables should be binary) is easily solved by simple code in a for loop:

Max:=-infinity:
for x__1 from 0 to 1 do
for x__2 from 0 to 1 do
for x__3 from 0 to 1 do
for x__4 from 0 to 1 do
for x__5 from 0 to 1 do
if 3*x__1+5*x__2+6*x__3+2*x__4+x__5 <= 10 then M:=3*x__1+14*x__2+18*x__3+6*x__4+2*x__5;
if M>Max then S:=[[x__1,x__2,x__3,x__4,x__5], M]; Max:=M fi;fi;
od: od: od: od: od:
S;

                                                  [[0, 0, 1, 1, 1], 26]

Here's another way to solve the same problem (more compact method suitable for more variables):

P:=combinat:-permute([0$5,1$5], 5):
Z:=(x__1,x__2,x__3,x__4,x__5) -> 3*x__1+14*x__2+18*x__3+6*x__4+2*x__5:
f:=(x__1,x__2,x__3,x__4,x__5) -> 3*x__1+5*x__2+6*x__3+2*x__4+x__5:
[seq(`if`(f(p[])<=10, [p, Z(p[])], NULL), p=P)]:
sort(%, (a,b)->a[-1]>b[-1])[1];

                                                   [[0, 0, 1, 1, 1], 26]

The third workaround is Mariusz Iwaniuk's way from here to increase  Digits :

restart;
Digits:=20:
Optimization:-LPSolve(3*x__1+14*x__2+18*x__3+6*x__4+2*x__5, {3*x__1+5*x__2+6*x__3+2*x__4+x__5 <= 10}, x__1 = 0 .. 1, x__2 = 0 .. 1, x__3 = 0 .. 1, x__4 = 0 .. 1, x__5 = 0 .. 1, assume=binary, maximize);

                           [26, [x__1 = 0, x__2 = 0, x__3 = 1, x__4 = 1, x__5 = 1]]

Edit.

Maybe you want it:

A:=Matrix(2, 2, {(1, 1) = (1/6)*sqrt(3)+(1/2)*I, (1, 2) = (1/6)*sqrt(3)-(1/2)*I, (2, 1) = (1/6)*sqrt(3)-(1/2)*I, (2, 2) = (1/6)*sqrt(3)+(1/2)*I});
A1:=(sqrt(3)/6)%*(Matrix(2,(i,j)->Re(A[i,j]))/(sqrt(3)/6))%+((1/2)%*(Matrix(2,(i,j)->A[i,j]-Re(A[i,j]))/(1/2)));   # Or
A2:=(sqrt(3)/6)%*(Matrix(2,(i,j)->Re(A[i,j]))/(sqrt(3)/6))%+((1/2)%*(Matrix(2,(i,j)->Im(A[i,j]))/(1/2))%*I);

If you want to calculate the values of an expression for different values of a variable, it is convenient to immediately specify it as a function, and not as an expression. Then in the future it will be very easy to calculate the value of this function symbolically or numerically.

Your example:
f := x->2*sin(x)-x^2/10:
f(5),  f(5.);
                                    2*sin(5)-5/2,  -4.417848549                    

Here are 3 ways to solve the problem by an example.

The problem: find the smallest root of equation  sin(100*x)+0.1*x=0  in the range  1 .. 2 . In this range there are several dozen roots.

The first way by  fsolve  command. Using  plot  command, we reduce the range to a single root:

restart;
eq:=sin(100*x)+0.1*x=0:
fsolve(sin(100*x)+0.1*x=0, x = 1 .. 2);
plot(lhs(eq), x = 1 .. 2);
plot(lhs(eq), x = 1 .. 1.1);
fsolve(lhs(eq), x = 1 .. 1.02);

The second way by  RootFinding:-Analytic  command, which provides all the real and complex roots in the range  re=1..2, im=-1..1  and we select the smallest real one:

[RootFinding:-Analytic(lhs(eq), re=1..2, im=-1..1)];
min(%);

The third way by  Student:-Calculus1:-Roots  command, which provides all the real roots in this range and we select the smallest one:

evalf(Student:-Calculus1:-Roots(lhs(eq), x = 1 .. 2));
min(%);