Use  plot3d  command instead of  plot  and unassigned variables  r  and  n :

plot3d(exp(-0.5*(r/sqrt((-r^2+1)/(n-2)))^2)/(sqrt((-r^2+1)/(n-2))*sqrt(2*Pi)), 'r' = 0.9 .. .1, 'n' = 4 .. 30, linestyle = solid, axes = normal, orientation = [-30, 70]);

In Maple 2015 and 2016  linestyle=solid  option is required:

plot3d(x^2+y^2, x = -1 .. 1, y = -1 .. 1, linestyle=solid);

Another way - not rewrite the system, but simply add extra equation in the original system:

eq1:=J*diff(theta(t),t,t)+b*diff(theta(t),t)=K*i(t):

eq2:=L*diff(i(t),t)+R*i(t)=V(t)-K*diff(theta(t),t):

eq3:=v(t)=diff(theta(t),t):

DCMotor:=[eq1,eq2,eq3];

plot([(x-1)/(x-2), [2, t, t = -5 .. 7], 1], x = -3 .. 7, -4 .. 6, color = [red, blue, blue], thickness = [2, 1, 1], linestyle = [1, 2, 2], discont, labels = [x, y]);

The same plotting for the second plot. Here the horizontal asymptote is  y=2  not y=1

compoly  is an alternative command for these simplifications:

compoly(numer(f), {x,y});

                                                 x^3-36*z^2, x = 5*x+7*y

That is, we have the identity

numer(f)=subs(%[2], %[1]);

                125*x^3+525*x^2*y+735*x*y^2+343*y^3-36*z^2 = (5*x+7*y)^3 - 36*z^2

The same for the denominator.

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