The following code is working:

restart; 
A := Matrix([[10, 2, 0, 0, 0], [1, 4, 2, 0, 0], [0, 1, 6, 2, 0], [0, 0, 2, 7, 3], [0, 0, 0, 2, 18]], datatype = float, storage = sparse);
b := Vector([11, 20, 34, 56, 77], datatype = float);
X := LinearAlgebra:-LinearSolve(A, b, method = SparseIterative);

Edit. The help says that the matrix can be nonsymmetric.

We see that this equation has 1 real and 2 complex roots:

[solve(5*x^3-x^2+x-1=0, explicit)];
evalf~(%);

Example:

restart;
P:=randpoly([x,y]);
C:=[coeffs(P,[x,y], 't')];
T:=[t];
TC:=zip(`[]`,T,C);
L:=[seq([[degree(p[1],x),degree(p[1],y)],p[2]], p=TC)];

I think this integral can only be calculated numerically when all parameters are given. 

Example:

restart;
Cap_Unit_Length := Int((epsilon_0*epsilon_r*d*sqrt(1+(T*Pi*S)^2))/(epsilon_r*S*theta*tan(theta/2)+2*t),theta=0..Pi);
evalf(eval(Cap_Unit_Length, [epsilon_0=1,epsilon_r=2,d=3,T=4,S=5,t=6]));

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]