@Axel Vogt Thank you for the simple way! :)

@Carl Love You and acer are totally right! :)
I refer you to my last answer here.

Thanks for the attention Carl. :)

@acer Actually, I am solving a boundary value problem by techniques of series expansions of sturm-liouville. First, I have to write the boundary conditions in the symbolic form, and then giving them numerical values for a specific boundary condition and specific number of tems included in the series solution. Therefore I can't write the data directly to the coefficient matrix "A" as someone in a FEM code usually do. For finding the coefficients of the series I have to solve a linear system. The linear system that I uploaded correspond to including 200 hundred terms in the series solution.

I was just asking myself that why the solution of linear system here takes lot more time than the FEM solution of the same problem. The answer to that was the numerical instability of my system. The presence of the terms such as cosh(z) and BesselI(0,r) in the series solutoin are the source of large numbers in the coefficients of the linear system as you have noticed but in the FEM analysis the coefficient matrix isn't ill-posed.

About the possiblity of the wrong generated system, I have verified my answer with a FEM solution of the same problem for 50 terms included in the seires solution. Also, I have checked it many times and I am sure the generated system is correct.

Now I am looking for a way to prevent the presence of large or small numbers in the coefficient matrix by some techniques. I hope that I proceed! :)

I really do appreciate your attention. As I told before, you are great! :)

@Kitonum Thank you very very much! You saved me from writting a procedure! :)

@acer Here I want to discuss about the idea of writting a procedure to develope the matrix form.

consider that you have a system of linear equations named as "EQ[i]" , i=1..N in a free from. Our goal is to produce the matrices "A" and "b". We know the unknown vector "x". For computing the components of "A", I will use the following:

EQ[i]:=op(1,EQ[i])-op(2,EQ[i])=0;
M:=Matrix(N,N,datatype=float[8]);
A[i,j]:=coeff(op(1,EQ[i]),x[j]);

How can I determine the components of the "b" matrix?

About "datatype=float[8]", is it used correctly here?

@acer
My system is not in a matrix form as you said. I will try to do this. I hope that this do the work with LinearAlgebra:-LinearSolve in few seconds. ;)

Have you looked at the file in the download link above? I will be so thankful if you do! :)
It just take a few seconds! ;)

@acer You are right about symbolic solution of a large system but my system is for a particular problem. My system doesn't contain any symbolic variable except the unknowns! How does the "solve" command try to solve such a system?

@acer By "purpose" I mean to solve the system as fast as possible! :)
I don't know that you are familiar with Finite Elment Softwares or not. In those sofwares the simulated problem reduces to the solution of large linear algebraic systems. Those sofwares like "ANSYS" or "ABAQUS" solve this large system in a minute! I don't know what they do! :D
But in "Maple" it takes hours! :)

@Alejandro Jakubi 

I couldn't understand that how the first way can be implemented here! Can you give me some hints? :)

About the second way, it's just like using a "subs" command since you can not pass the integrand to the "int" command directly! You just substitute some variables in a predefined formula! You can do it with just creating a procedure named "rBB" and not " `int/rBB` ". I mean that calling the "int" command in here is unnecessary! :)
See the file below:

Bessel_Int_Procedures.mw

What I am looking for is to be able to pass the integrand directly to the "int" command.

@acer I mean to do something that enable maple to compute these kinds of integrals automatically! :)
As you see in "Integral.mw" maple can not compute the anitderivative function. I want to help maple to compute this integral automatically. I was wondering that there should be a way!

@acer I have an integral that maple can not solve but I can solve it by hand. How can I add this to maple integration database?

f:=int(r^2*BesselJ(0,a*r)*BesselI(1,b*r),r)

Please see the file below.

Integral.mw

@acer You are great.
Did you do anything about the contour plot question mentioned here?

@Kitonum It is really better than imaginary exponentials. Thank you very much. Can you do something that exponentials be replaced by hyperbolic funtions "cosh(b*x)" and "sinh(b*x)"?

@Markiyan Hirnyk Continnue being glad! :p
All others have accepted that such a feature should be included and I don't know what's your problem glad guy! :D

@ecterrab That was so helpful. Thank you very much! :)