I think it can be done in different ways. Here's one way as an example:

A:=`01101001`:

`+`(convert(parse(A), base, 10)[]);

                                    4

 For your example:

mylist:=[4,6,2]:

z:=map(convert, mylist, binary):

map2(nprintf, "%08d", z):

map(t->`+`(convert(parse(t), base, 10)[]), %);

                                  [1,2,1]

I did not find any mistakes. Just your function decreases very rapidly. It can be clearly seen, if you reduce the range for x-axis:

plot((11), x=0..1, y=0..0.4, labels=["x(m)","PDF(x)"], color=red,thickness=3 );

Addition:

plot((11),x=0..1, y=0..7, labels=["x(m)","PDF(x)"], color=red,thickness=3 );

As OP wrote "the variables are u, v, w, alpha, beta and gamma"  first we define the vector  X(t)  that specifies the order of unknown functions, and then subsequently find matrices  M, C, K  and vector  F :

V:=Vector(6, {(1) = (2*R^2*(diff(w(t), t))*Pi*Omega*h*rho+R^2*(diff(u(t), t, t))*Pi*h*rho-R^2*u(t)*Pi*Omega^2*h*rho = 0), (2) = (R^2*(diff(v(t), t, t))*Pi*h*rho = 0), (3) = (-2*R^2*(diff(u(t), t))*Pi*Omega*h*rho+R^2*(diff(w(t), t, t))*Pi*h*rho-R^2*w(t)*Pi*Omega^2*h*rho = 0), (4) = ((1/4)*R^4*Pi*(diff(alpha(t), t, t))*h*rho+(1/12)*R^2*Pi*(diff(alpha(t), t, t))*h^3*rho+(1/6)*R^2*Pi*(diff(gamma(t), t))*Omega*h^3*rho-(1/12)*R^2*Pi*alpha(t)*Omega^2*h^3*rho = 0), (5) = ((1/2)*R^4*Pi*(diff(beta(t), t, t))*h*rho-(1/2)*R^4*Pi*beta(t)*Omega^2*h*rho = 0), (6) = ((1/4)*R^4*Pi*(diff(gamma(t), t, t))*h*rho+(1/12)*R^2*Pi*(diff(gamma(t), t, t))*h^3*rho-(1/6)*R^2*Pi*(diff(alpha(t), t))*Omega*h^3*rho-(1/12)*R^2*Pi*gamma(t)*Omega^2*h^3*rho = 0)});

X(t):=<u(t), v(t), w(t), alpha(t), beta(t), gamma(t)>;

M:=Matrix(6, (i,j)->coeff(lhs(V[i]), diff(X(t)[j],t,t)));

C:=Matrix(6, (i,j)->coeff(expand~(lhs~(V) - M.map(diff, X(t),t,t))[i], diff(X(t)[j],t)));

K:=Matrix(6, (i,j)->coeff(expand~(lhs~(V) - M.map(diff, X(t),t,t) - C.map(diff, X(t),t))[i], X(t)[j]));

F:=-expand~(lhs~(V) - M.map(diff, X(t),t,t) - C.map(diff, X(t),t) - K.X(t));

Addition: if the order of unknown functions is specified (vector X(t)), we can easily, using a for loop, write a procedure that generalizes the solution to any linear system.

System.mws 

Should be:

ContoursWithLabels(PPP, -1 .. 200, -1 .. 110, [thickness = 3]);

A version with coloring and with selected contours :

ContoursWithLabels(PPP, -1 .. 200, -1 .. 110, {seq(10^(-5) .. 2.2*10^(-4), 3*10^(-5))}, [color = black, thickness = 2, view = [-1 .. 200, -1 .. 110]], Coloring = [colorstyle = HUE, colorscheme = ["LightYellow", "Orange", "Red"], style = surface]);

 Edited.

Your equation has a unique root in any range  -Pi/20+Pi*n/10<=sigma<=Pi/20+Pi*n/10 , n is an integer.

An example for n=1:

fsolve(tan(10*sigma)+tanh(10*sigma) , sigma=Pi/20..3*Pi/20);

                                             0.2365020372

The method is based on direct calculation of  1/p , and can be applied to any fraction  1/p   in which  p is prime and p<>2 and p<>5 :

myproc1:=proc(p::prime)

local n;

n:=numtheory[order](10, p);  # Lenght of the period

evalf(1+1/p, n+2):

convert(%, string):

n, %[3..n+2];  # Lenght of the period and the period itself

end proc:

Examples of use:

myproc1(7);

myproc1(17);

myproc1(37);

myproc1(97);

In fact, we want to create a vector function of a scalar argument, from which it is easy to extract individual  "coordinate" functions.  

An example  (I use a list instead of Array or Vector):

A:=[sin(x), x^2, x^3-2];

A:=unapply(A, x);

A(m)[1], A(m)[3];

a:=[1,2,3,4,5]:

convert(cat(op(a)), string);

                                               "12345"

applyrule(exp(a::anything + lambda::realcons*b::symbol) = exp(a)*exp(lambda*b), F);

If the series is absolutely convergent, then it can be summed by a blocks, for example

evalf(sum(sum(1/(i^4+j^4), i=0..infinity), j=1..infinity));

                                        1.844054925+0.*I

A := plot(x^2, x = -2 .. 2, color = red, thickness = 2, axis[2] = [color = white], labels=[" "," "]):

B := plots[textplot]([[1.9, -0.1, "x values"], [0, 4, "y values"]], align = [below, left], font = [TIMES, ROMAN, 14]):

plots[display](A, B);

Edited.

The procedure  IntersectionOfSegments  solves the problem for any segments on the plane. If the segments do not intersect, then the procedure returns NULL. If segments intersect at the single point, the procedure returns this point (as a list). If the intersection of segments is a segment, then the procedure returns this segment (as a listlist).

IntersectionOfSegments:=proc(S1::listlist,S2::listlist)

local A, B, C, D, f, g, T1, T2, S11, S21, T;

uses geometry;

point(A, S1[1]); point(B, S1[2]); point(C, S2[1]); point(D, S2[2]);

f:=unapply(lhs(Equation(line(AB, [A,B]),[x,y])),x,y);

g:=unapply(lhs(Equation(line(CD, [C,D]),[x,y])),x,y);

if is(f(S2[1][])*f(S2[2][])<=0) and is(g(S1[1][])*g(S1[2][])<=0) then  if is(f(S2[1][])<>0) or is(f(S2[2][])<>0)  then return rhs~(convert(solve({f(x,y), g(x,y)}),list))  elif

is(S1[1,1]<>S1[2,1]) then S11:=sort(S1,(a,b)->is(a[1]<=b[1])); S21:=sort(S2,(a,b)->is(a[1]<=b[1])); if is(S11[2,1]>=S21[1,1]) and is(S11[1,1]<=S21[2,1]) then T1:=sort([S11[1],S21[1]],(a,b)->is(a[1]<=b[1])); T2:=sort([S11[2],S21[2]],(a,b)->is(a[1]<=b[1]));

T:=[T1[2],T2[1]];

return `if`(T1[2]<>T2[1],sort([T1[2],T2[1]],(a,b)->is(a[1]<=b[1])), T1[2]) fi else

S11:=sort(S1,(a,b)->is(a[2]<=b[2])); S21:=sort(S2,(a,b)->is(a[2]<=b[2]));

if is(S11[2,2]>=S21[1,2]) and is(S11[1,2]<=S21[2,2]) then T1:=sort([S11[1],S21[1]],(a,b)->is(a[2]<=b[2])); T2:=sort([S11[2],S21[2]],(a,b)->is(a[2]<=b[2]));

T:=[T1[2],T2[1]];

`if`(is(T1[2]<>T2[1]),sort([T1[2],T2[1]],(a,b)->is(a[2]<=b[2])), T1[2])

fi;  fi; fi; 

end proc:

Examples of use:

IntersectionOfSegments([[0,1],[1,0]], [[0,-1],[1,2]]);

IntersectionOfSegments([[0,1],[1,0]], [[1/2,1/2],[2,-1]]);

IntersectionOfSegments([[0,1],[1,0]], [[-2,3],[-1,2]]);

IntersectionOfSegments([[0,0],[0,2]], [[0,1],[0,3]]);

IntersectionOfSegments([[0,0],[0,2]], [[0,5/2],[0,3]]);

IntersectionOfSegments([[0,0],[0,2]], [[0,2],[0,3]]);

 Segments.mws

I have plotted your system  for several values of  r, from which we can see:
1) There are always 2 positive solutions (2 points in the ranges omega=0..1, a=0..1).
2) These two points are located far enough from one another.

Further, using  fsolve  command, we find all the solutions for  20 values  r=0.1 .. 2

We have plotted only 4 pairs of points as for  r>=1 the points practically merge.

Download Syst1.mw

Your integrand is only defined on a finite set of points. To integrate, it is necessary  to continue on a certain range. For example on the range t=1..10  you can use splines:

x:=[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]:

rho:= [0.045, 0.0459, 0.0564, 0.05689, 0.06015, 0.06235, 0.0654, 0.0687, 0.07012, 0.07251]:

eta:= [1.15, 1.256, 1.56, 1.85, 1.86, 2.01, 2.35, 2.56, 2.86, 2.901]:

f:=[seq(rho[i]*x[i]^4*eta[i]^2, i=1..10)]:

F:=unapply(CurveFitting[Spline](x,f,t), t);

plot(F,1..10, color=red, thickness=2);

int(F(t), t=1..10);

 Addition. For some reason, Maple does not plot for t>8, although at certain points value of the function is calculated correctly:

F(9), F(10);

f[9], f[10];

                               3763.08485500000, 6102.29730500000

                                        3763.084855, 6102.297305

A bug?

p  is a polynomial not necessarily a cubic polynomial.

CriticalPoints:=p->solve(diff(p, x), x):

Example of use:

CriticalPoints(x^2*(x-1));

                                                     0, 2/3