If you're going to call the function z 800 times, it looks very inefficient to call CDF twice in each of those calls to z, rather than just making two of those calls in all.  Thus you might try something like this.

> FN:= unapply(CDF(Normal(0,3),x), x);
   FG:= unapply(CDF(Gumbel(0,1), y), y);
   z:= (X,Y) -> min(FN(X), FG(Y));

 Also, did you really mean to use Gumbel(0,1) in one place and Gumbel(0,2) in the other?

You can't "get rid of it", but you can evaluate it numerically, like most other special functions. 

> evalf(Psi(1/2 + 1/6*sqrt(6)), 50);

−0.73917729505204124955610221184504282302704283831579

You seem to have a "cloud" of 1000 points in 3 dimensions.  A 2-d contourplot is the plot of some function f(x,y) of two variables.  What do you want this function to represent?

The symbolic differential equation solver seems to have trouble with this one, so you could try the numeric solver.

> sop:= dsolve({u,ins}, {A(t),B(t),F(t),M(t),K(t)}, numeric);

> plots[odeplot](sop, [t, A(t)], t = 0 .. 25);

I didn't go up to 1000, because the solution "blows up" long before that.

In fact, what happens here is that the coefficient matrix for this linear system has two positive eigenvalues, 
approximately 0.791601968171963266 and 0.135156910763671200, which will dominate the solution for large t.

You could also try using LinearAlgebra[MatrixExponential]: if you write the system in vector form as

D(V)(t) = Q.V(t) + b, then the solution is

V(t) = MatrixExponential(Q*t) . (V(0) + Q^(-1).b) - Q^(-1).b

Thus:

> with(LinearAlgebra):  Digits:= 14;
    Q,b := GenerateMatrix(subs(A(t)=a,F(t)=f,B(t)=b,K(t)=k,M(t)=m,map(rhs,[u])), [a,b,f,m,k]);
    U:= MatrixExponential(Q,t);
    V0:= subs(ins,<A(0),B(0),F(0),M(0),K(0)>);
    S:= U.(V0 - Q^(-1).b) - Q^(-1).b;

> plot(S[1], t = 0 .. 25);

evalm and &* are for the old-style matrix type, not for Matrix.  Why are you using them?

It is rather strange that the help page uses Matrix rather than matrix. 

For the domain, think: what values of x would make this function undefined?  Typical problems: division by 0, square root or logarithm of a negative number (assuming complex numbers are not allowed).

Range can be more difficult. Plotting the function is often useful.

The GF package is rather old, and doesn't seem to interact well with linear algebra.  Can you avoid using it?

For example, to work in a field of order 2^4:

> P:= x^4 + x + 1;
Check that this is irreducible mod 2:

> Irreduc(P) mod 2;

       true

> alias(alpha = RootOf(P, x));

> with(LinearAlgebra):
   M:= << 1 + alpha + alpha^2, 1 + alpha + alpha^3> | <0, 1 + alpha >>;
   V:= < 1 + alpha^2, alpha + alpha^3>;

Then to multiply these:

> simplify(M . V) mod 2;

For example, to start you off:

> with(plots): with(Statistics):
  sigma:= 1: mu:= 1:
  randomize():

Standard Brownian motion W(t) may be approximated using a random walk where each step over time interval h is a normal random variable with mean 0 and standard deviation sqrt(h).

>  h:= 0.001: 
   Y:= RandomVariable(Normal(0,sqrt(h)));
   for k from 1 to 10 do
     L:= Sample(Y,1000);
     B[0]:= 0; 
     for i from 1 to 1000 do B[i]:= B[i-1]+L[i] end do:
     pts[k]:= <seq(<h*j, exp((mu-sigma^2/2)*h*j + sigma*B[j])>,j=0..1000)>:
   end do: 
   display([seq(pointplot(pts[k], style=line),k=1..10), 
      plot(exp(mu*t), t=0..1)]);

Make a list of the real and imaginary parts, and use pointplot.

In general, to apply some function to all elements of a Matrix, list or set, use map.

> H:= <1 | x | x^2 | x^3>;
  B:= <1,0,0,0|0,1,0,0|-3/L^2, 2/L, 3/L^2, -1/L| 2/L^3, 1/L^2, -2/L^3, 1/L^2>;
  map(int, H.B, x = 0 .. L);

[L, 1/2*L^2, -3/L+2*L-1/4*L^3, 2/L^2-1/6+1/4*L^2]

Make it into a function of f.  If you already have the equation shown (I'll call it eq), then

> res:= unapply(rhs(eq), f);

You will probably want to use a logarithmic scale for the frequency, otherwise all the action will be in one tiny corner of the graph.

> plots[semilogplot](res(f), f = 1 .. 10^9);

Of course the other parameters in your equation, a, La, rho, mu[r] and mu[0], will need to have numeric values for this to work.

To start you out: you might try this:

> with(plots):
   P1:= plot3d(1, phi=0..2*Pi,theta=0.. arccos(-1/3), coords=spherical, style=wireframe);
   P2:= plot3d(1, phi=0..2*Pi,theta = arccos(-1/3) .. Pi, coords = spherical, style=patchnogrid):
   display([P1,P2]);
 

 For the disks, you might find cylinder in the plottools package useful.  Or perhaps VolumeOfRevolution
in Student[Calculus1].

That should be Fit rather than fit.  Maple is case-sensitive.

It seems that a bug in the 2D Input parser is preventing this from working.  But why are you trying to do this with "use"?  It's very simple to
define nested procedures without "use".  Just assign the nested procedure as a local variable.  For example:

plotHCDF= proc(x, a, b, itype, x1, x2)
local barArea, xmin, xmax, p, H1, H2, m, s, SL, H, A;
use Statistics, ... in
  barArea:= proc(polygon::listlist, $):: nonnegative;
      ...
    end proc;
  ...
end proc;