You can easily do this by removing the term  O(x^6)  and adding the required term (here  xi  is between  0  and  x ):

restart;
convert(series(cos(x),x=0, 5), polynom)+(D@@6)(cos)(xi)*x^6/6!;

See the general formula in  https://en.wikipedia.org/wiki/Taylor%27s_theorem

Edit.

I think that different forms of answers are related to different methods used in the calculation. Use method=FTOC option for simplicity in this example:

int(2*Pi*x*sqrt(1+x^(-4/3)/9), x=0..1, method=FTOC);

FTOC (abbreviation) = fundamental theorem of calculus

Since the interval for the unknown  h  (in which we are looking for a solution) has already been indicated  h>=10 and h<=20, everything can be done in 1 step:

restart;
f := x -> sqrt(-x^2+20*x):
solve(2*Pi*int(f(x)*sqrt(1+diff(f(x),x)^2),x=0..h)=1005) assuming h>=10,h<=20;

                                         201/(4*Pi)

We can easily find all matrices that commute with A by simply reducing the problem to solving a system of 4 equations with 4 unknowns:

restart;
A := Matrix(2, 2, [[3, 1], [-3, 2]]):
C := Matrix(2, 2, [[x, y], [u, v]]);
A.C=~C.A;
Sys:=convert(%,set);
solve(Sys);
C:=eval(C, %); # The general formula for C,where  x and y are any numbers
E:=eval(C, [x=1,y=2]); # The specific example
E.A=A.E;  # Check

The new procedure  NewLegend1  allows you to place the opaque rectangle anywhere on the plot. Please note that in the list  G , you must specify a spacing for the variable  y  also as for  x .

restart;
NewLegend1:=proc(f::algebraic,g::algebraic,R::list,G::list)
local A1, A2, A3, A4, B, a, b, l1, l2, T;
uses plottools, plots;
A1:=implicitplot([y=f,y=g], subsop(1=(x=op([2,1],G[1])..R[1,1]),2=(y=op([2,1],G[2])..R[1,2]),G)[], labels=["x","y"]);
A2:=implicitplot([y=f,y=g], subsop(1=(x=R[1,1]..op([2,2],G[1])),2=(y=op([2,1],G[2])..R[2,2]),G)[], labels=["x","y"]);
A3:=implicitplot([y=f,y=g], subsop(1=(x=R[2,1]..op([2,2],G[1])),2=(y=R[2,2]..op([2,2],G[2])),G)[], labels=["x","y"]);
A4:=implicitplot([y=f,y=g], subsop(1=(x=op([2,1],G[1])..R[2,1]),2=(y=R[1,2]..op([2,2],G[2])),G)[], labels=["x","y"]);
B:=rectangle(R[1],R[2], style=line);
a:=R[2,1]-R[1,1]; b:=R[1,2]-R[2,2];
l1:=line(R[1]+[0.1*a,-0.3*b],R[1]+[0.4*a,-0.3*b],color=rhs(G[3])[1],thickness=rhs(G[4]));
l2:=line(R[1]+[0.1*a,-0.7*b],R[1]+[0.4*a,-0.7*b],color=rhs(G[3])[2],thickness=rhs(G[4]));
T:=textplot([[(R[1]+[0.75*a,-0.3*b])[],y=f],[(R[1]+[0.75*a,-0.7*b])[],y=g]], font=[times,bold,14]);
display(A1,A2,A3,A4,`if`(nops(R)=2,B,NULL),l1,l2,T);
end proc:

Example of use:

NewLegend1(sin(x),cos(x), [[0.2,0.9],[1.4,0.3]], [x=0..2*Pi, y=-1..1, color=[red,blue], thickness=2, size=[900,400], scaling=constrained, axes=box]);

The  NumberTheory:-Divisors  command returns the set of all divisors of an integer. Procedure below solves your problem:

restart;
a:=proc(n)
local S, N, i, j, k;
uses NumberTheory;
k:=0;
S:=NumberTheory:-Divisors(n); N:=nops(S);
for i from 1 to N-2 do
if ormap(j->irem(S[j],S[i])=0, [$i+1..N-1]) then k:=k+1; fi;
od:
k;
end proc:

seq(a(n), n=1..24);

                               0, 0, 0, 1, 0, 1, 0, 2, 1, 1, 0, 3, 0, 1, 1, 3, 0, 3, 0, 3, 1, 1, 0, 5

Edit. The procedure  a  can be written much shorter:

a:=n->nops(NumberTheory:-Divisors(n))-nops(NumberTheory:-PrimeFactors(n)) - 1:

The third way: the procedure  a  can be written without calling the  NumberTheory  package, only using the kernel command  ifactors :

a:=proc(n)
local L, N;
L:=ifactors(n)[2];
N:=nops(L);
mul(L[i,2]+1, i=1..N)-N-1;
end proc:

Below is a possible variant of the desired procedure. The list  R  defines the position of the rectangle with legends. If the third element of the list  R  is some symbol, for example 'off', then the rectangle is not drawn.

restart;
NewLegend:=proc(f::algebraic,g::algebraic,R::list,G::list)
local A, B, a, b, l1, l2, T;
uses plottools, plots;
A:=plot([f,g], G[], labels=["x","y"]);
B:=rectangle(R[1],R[2], style=line);
a:=R[2,1]-R[1,1]; b:=R[1,2]-R[2,2];
l1:=line(R[1]+[0.1*a,-0.3*b],R[1]+[0.4*a,-0.3*b],color=rhs(G[2])[1],thickness=rhs(G[3]));
l2:=line(R[1]+[0.1*a,-0.7*b],R[1]+[0.4*a,-0.7*b],color=rhs(G[2])[2],thickness=rhs(G[3]));
T:=textplot([[(R[1]+[0.75*a,-0.3*b])[],y=f],[(R[1]+[0.75*a,-0.7*b])[],y=g]]);
display(A,`if`(nops(R)=2,B,NULL),l1,l2,T);
end proc:

Examples of use:

NewLegend(sin(x),cos(x),[[3.5,0.9],[4.3,0.3]],[x=0..2*Pi,color=[red,blue],thickness=2,size=[900,400],scaling=constrained, axes=box]);

NewLegend(exp(x),x,[[2.5,1.8],[3.6,1],'off'],[x=-2..5,color=[red,blue],thickness=2,size=[500,500],scaling=constrained, view=[-2..5,-2..5],axes=normal]);

Edit. 1. I don't understand why you color the rectangle and place it in the top left corner every time. If you want the procedure to be universal, you cannot specify in advance the place where the rectangle should be located. It may turn out that just in this place there will be something interesting on the plot, and the rectangle will block it. Also, with your way, you won't be able to use the  view  option which is often useful (see your second example with  exp(x) ).
2. The procedure is written for the case of two graphs, but it is easy to modify it to suit any number of graphs, for example, if you have only one function or 3 functions, etc.
3. Please note that all the options required for plotting are included in the list  G . I understand your desire to reduce the number of formal parameters, but the degree of automation should be reasonable. The main thing in my opinion is the best end result.

To extract single entries or entire rows or columns use square brackets instead of parantheses:

restart;
A:=Matrix([[1,2],[3,4]]);
A[1,..];
whattype(%);
A[..,1];
whattype(%);

For rows  A[1,..] = A[1]  and so on.

It looks like

restart;
f:=(x,t)->sqrt(2)*exp(-0.6*(0.5*x+0.01*t)^2):
plot([f(x,0),f(x,20),f(x,40)], x=-10..10, color=[yellow,black,blue], thickness=3, size=[900,500]);

What 2 lines did you mean? In the plot from your worksheet, I see only one straight line and 2 coordinate axes.
General advice: to find the intersection point of two lines  F(x,y) = 0  and  G(x,y) = 0 , solve the system

solve({ F(x,y) = 0,  G(x,y) = 0}, {x,y});

Edit. Sorry, I didn't notice your second plot right away. So do as Carl wrote or just

restart;
solve({y=5*x-1, y=-x+7});

restart; 
P := 2^n-1: 
select(isprime, [seq(eval(P, n = k), k = 1 .. 100)]);

  [3, 7, 31, 127, 8191, 131071, 524287, 2147483647, 2305843009213693951, 618970019642690137449562111]

Edit. If you do not need prime numbers of the form  P  themselves, but exponents  n , then below is another code:

restart; 
P := 2^n-1: 
select(p->isprime(eval(P, n=p)), [$1..100]);

                 [2, 3, 5, 7, 13, 17, 19, 31, 61, 89]

restart;
a:=3: R:=7: z:=sqrt(R^2-(r-a+R)^2):
plot3d([[r,phi,z],[r,phi,-z]], phi=0..2*Pi, r=2..3, coords=cylindrical, scaling=constrained);

The formula  sqrt(R^2-(r-a+R)^2)  was obtained by elimination  phi  from the system  { z=R*cos(phi), r=a-R*(1-sin(phi)) } 

Edit. I added the option  scaling=constrained so that you can watch the form change depending on the parameters change.

In my opinion the simplest way is to use a triple for-loop:

restart;
C := y^2*z + y*z^2 = x^2:
for x in {0,1} do
for y in {0,1} do
for z in {0,1} do
if C then print(['x','y','z']=[x,y,z]) fi;
od: od: od:

To plot a 3-variable equation, you must use  plots:-implicitplot3d  command. To get the best quality we can first express  z  through  x  and  y and then plot:

restart;
eq:= x^2-y^2*z-y*z^2;
plots:-implicitplot3d(eq,x=0..1,y=0..1,z=-7..7, style=surface, grid=[40,40,40]);
solve(eq, z);
plot3d([%], x=0..1,y=0..1, style=surface);

Edit. A more interesting picture is obtained if we increase the ranges for  x  and  y :  x=-1..1, y=-1..1. It is worth noting one drawback of  plot3d  command in this example: Maple does not draw a surface near z-axis. This is due to the fact that for small values of  y , the value of  z  tends to infinity. The use of  plots:-implicitplot3d  command does not have this drawback. We get good quality when used together  plot3d  and  plots:-implicitplot3d :

restart;
eq:= x^2-y^2*z-y*z^2;
A:=plots:-implicitplot3d(eq,x=-1..1,y=-1..0.01,z=-5..5, style=surface, grid=[50,100,50]):
Sol:=solve(eq, z);
B:=plot3d([Sol], x=-1..1,y=-1..1, view=-5..5, style=surface, axes=normal, grid=[500,500]):
plots:-display(A,B, view=[-1..1,-1..1.2,-5..5]);

Let's help Maple a little. If  t^2 - 4 = s^2 - 4  and  t<>s  then  s = -t . Below we get that there will be infinitely many double points:

restart;
x:=t->sin(t): y:=t->t^2-4:
solve({sin(t)=sin(s),s=-t}, AllSolutions);
eval(%, _Z1=k);
[x,y]=eval([x(t),y(t)], %) assuming k::posint;  # All the double points
# k = 1, 2, 3, ...

                               [x, y] = [0, Pi^2*k^2-4]

The numerical solution:

restart;
x:=t->sin(t): y:=t->t^2-4:
Sys:=[x(t)=x(s),y(t)=y(s)]:
plots:-implicitplot(Sys, t=-10..10,s=-10..10, color=[red,blue], gridrefine=3);
fsolve(Sys, {t=-4..-1,s=1..4});
fnormal(eval([x(t),y(t)], %));

                                    [0., 5.869604404]

Because there are many solutions, then the plot is needed to isolate the root.