(1) use " instead of ', e.g

   RETURN ("Input is NOT a Lie algebra (',i,i,k,')=",A[i,i,k], "is not zero"); 

(2)

for l from 1 by 1 to n do

(you have the number "one" instead of letter l )

(edited)

flist := [0.1,0.2,2,3,3.5,6]:
xlist := [1,2,3,4,5,6]:
add( (xlist[i]-xlist[i-1])*flist[i]*(xlist[i]^2-1), i=2..6 );
                             355.6

or

add( (xlist[i+1]-xlist[i])*flist[i]*(xlist[i]^2-1), i=1..5 );
                             145.6

a2:=(18*d^2*R^2*r-4*d1^3*r+5*R^6-z1*r^3*p-14*r1^3*r*rc+k1*d^2*R^2*r);
select(hastype,[op(a2)],negint); a3 := add(i,i=%);
a4 := a2-a3;

Do you mean that A contains the first en B the second argument of BesselY?
In that case (I changed B a little):

A := Matrix([[-13, -10], [21, 16]]): B := Matrix([[3, 1], [21, 16]]):
BesselY~(A,B);

You forgot some multiplication signs inthe second equation:

restart;
eq1:=y=86:
eq2:=y=-0.0000054527*x^3+0.010903836*x^2+0.0714244709*x+74.18816:
sol:=solve({eq1,eq2},{x,y});
    {x = 30.00404568, y = 86.}, {x = 2005.705515, y = 86.},  {x = -35.99639135, y = 86.}

evalf( Int( 2.91*x*((1/(1+1.38*x^4)))^0.431 - 3.459*x^5/((1/(1+1.38*x^4))^0.569*(1+1.38*x^4)^2), x=0..1 ) );
                          1.001485791

ShowCols := proc(A)
  print( seq(LinearAlgebra:-Column(A,i), i=1..LinearAlgebra:-ColumnDimension(A)) )
end proc:

with(LinearAlgebra):
A := RandomMatrix(5,7);
ShowCols(A);

plots:-pointplot([seq([x,0],x=[-sqrt(16),-2,-3/4,0,exp(1),Pi])],view=[-4..4,0..0.01],
  symbolsize=20,symbol=solidcircle,colour=blue, 
  tickmarks=[[-4.0=typeset(sqrt(`16`)), -2.=-2, -.75=typeset(-3/4),0.=0,
  evalf(exp(1))=typeset(exp(1)),evalf(Pi)=typeset(Pi)],[]], scaling=constrained);

You must ensure that the expression contains plus-signs, otherwise there is only one term. What about:

f := (b(t)*diff(a(t),t) + a(t)*b(t))*k ;
if type( expand(f),`+`) then nops(expand(f)) else 1 end if; #corrected

Solve Im(SS)=0, use the avoid option to get all solutions:

z1 := fsolve( Im(eval( SS, a=0.3 ))=0, {_Z2} );
                           {_Z2 = 0.}
z2 :=fsolve( Im(eval( SS, a=0.3 ))=0, {_Z2} , avoid=z1 );
                          {_Z2 = -1.}
z3 :=fsolve( Im(eval( SS, a=0.3 ))=0, {_Z2} , avoid=( z1 union z2) );
      fsolve(-Im(LambertW(_Z2, -0.7 exp(-1))) = 0, {_Z2},  avoid = {_Z2 = -1., _Z2 = 0.})

No further solutions.

L := [1,2,3,7,6,5,4]:
m,pos := ListTools:-FindMaximalElement(L,position);

If you want to assume that a and b have the same sign, you can substiture b = k*a, and assume that k>0:

eq :=(4*a^3*b)^(1/2)/(-(a/(4*b))^(1/2))+(4*a^3*b*(4*b/a))^(1/2) = 0:
eq1 := eval(eq, b=k*a ):
simplify( eq1 ) assuming k>0 ;

You could use MathType to convert LaTeX code to Word.

No chance to find exact solutions as af function if f1. So use fsolve for each desired value of f1:

ABC := t -> fsolve( eval({A,B,C},f1=t), {a,b,c} ):
ABC(0.1);

And to get a as a function of t:,

AA := t -> subs( ABC(t),a );

but the solutions seen rather chaotic (I did not check your formulas)

The most elegant way to produce this kind of matrices is to use an indexing function:

A := M -> Matrix(2*M+1, (i,j) -> if (i=M+1 and j=M+1) then 1
                                 elif j=M+1 then 1/(i-1-M)/Pi 
                                 elif i=M+1 then 1/Pi
                                 elif (i+j)=2*M+2 then 1/abs(i-1-M)/Pi 
                                 else 0 
                                 end if):
A(3);