To generate a sample of size 99, for example, do

X:= Statistics:-RandomVariable(ProbabilityTable([.2, .5, .1, .1, .1])):
S:= Statistics:-Sample(X, 99):
trunc~([seq(S)])

If the support is not an initial segment of positive integers, then use EmpiricalDistribution instead of ProbabilityTable.

There's a factor of 1/sqrt(x) at the start of the expression. Thus, you need to use 

limit(asa, x= 0)

rather than

subs(x= 0, asa)

Do

remove(type, ListTools:-Flatten(w), {constant, linear})

Or, if you want to preserve the sublist structure (no matter how complicated), do

evalindets[2](w, list, type, remove, {constant, linear});

In either case, the resulting list of nonlinear terms can be counted by

nops(ListTools:-Flatten(%))

Maple has no problem doing the animation. The problem was you not stating the problem clearly. 

My animation is similar to Acer's, but I used a tangent transformation to smooth the rotation. 

restart:
Cons:= {4*x1+x2 <= 12, x1-x2 >= 2, x1 >= 0, x2 >= 0}:
Obj:= c1*x1 + 2*x2: 
LP:= proc(C1)
option remember;
    if C1::numeric then
         Optimization:-LPSolve(eval(Obj, c1=  C1), Cons, maximize)[1]
    else
        'procname'(C1)
    fi
end proc
:
plots:-animate(
    plot, 
    [
        (LP(tan(c1))-tan(c1)*x1)/2, x1= 0..4, x2= -1..2, thickness= 3, color= red,
        title= 'sprintf'("c1 = %4.2f, Z = %4.2f", tan(c1), LP(tan(c1)))
    ],
    c1= -Pi/2..Pi/2, frames= 300, paraminfo= false,
    background= plots:-inequal(Cons, x1= 0..4, x2= -1..2, nolines)
);

There's no way that you're going to get an interesting result for this, because Maple doesn't even know how to differentiate KummerU with respect to its first parameter. But if you just want a result, any result, instead of an error message, you could do

series(' 'KummerU' '(p, 3/2, t), p);

That's two set of single quotes, not one set of double quotes.

Your underlying plot command is nonsense because it specifies a plot with respect to two variables, x1 and x2. There can only be one such variable.

Your trouble has nothing to do with the background option.

You can't specify both the value of the step size h and the desired accuracy at the same time because they depend on each other.

You'll also need to specify either the number of steps or the final x-value. Usually n is the number of steps, so are you sure that you're interpretting the instructions correctly?

Your dx and point_list don't make sense because they're loop invariants: Their value doesn't change in the loop.

I do this:

solve({sec(x)=sqrt(2), x>=3*Pi/2, x<=2*Pi}, allsolutions, explicit);

This'll also work with RealDomain:-solve. IMO, it's usually better to put inequalities inside the solve rather than in an assuming clause. (And the help ?RealDomain says that assumptions don't work with that package.)

This usage of explicit is undocumented, but I've been using it since Maple 16 (and it may have existed before that). I just discovered it by experimentation. Both the keywords allsolutions and explicit are required to get the specific results that satisfy the inequalities. Taking the help page ?solve,details at face value, it makes no sense to me that these keywords would work together for a transcendental equation, but the above results show that they obviously do. 

I would avoid using RealDomain unless you're sure that you need it. Wanting only real solutions is not a good reason to use it.

You can provide fsolve with ranges within which to search, if you know some:

fsolve({r1, r2}, {xa1= 0..1, xb1= 0..1});
      {xa1 = 0.5048762298, xb1 = 0.06850883394}

This solution is returned immediately.

Any system of recurrence relations of finite order with initial conditions (regardless of linearity, homogeneity, or autonomy) is equivalent to an ODE IVP system where the first derivatives are interpretted as the discrete derivative diff(x(n), n) = x(n+1) - x(n), and this is extended to higher-order derivatives in the natural way: (D@@k)(x)(n) = (D@@(k-1))(x)(n+1) - (D@@(k-1))(x)(n). The new equations are formally called difference equations rather than differential equations. If that ODE IVP is passed to dsolve(..., numeric, ...), and solved with Euler's method with stepsize = 1, then the solutions are identical to the solutions of the recurrence relations. These solutions can be plotted with plots:-odeplot, which makes this a very convenient method.

Here's a utility to convert a system of recurrence relations to its ODE IVP equivalent:

RecToDsys:= proc(
    Sys::set(algebraic=algebraic), #recurrence relations expressed with indexed names
    V::list(name), #dependent variables
    n::name, #independent variable
    ICs::set(indexed(integer)=complexcons) #initial conditions
)
local Z,K,k,j;
    (Z,K):= (min,max)(op~(lhs~(ICs)));
    evalindets(
        eval(Sys, n= n-min(eval(op~(indets(Sys, specindex(V))), n= 0))),
        specindex(V), 
        t-> 
            local f:= op(0,t), k, K:= eval(op(t), n= 0);
            add(binomial(K,k)*diff(f(n), [n$k]), k= 0..K)
    ) union 
    eval(
       `union`(
            seq(
                {(
                    (D@@k)~(V)(Z)=~ 
                    `~`[`+`](seq((-1)^(k-j)*binomial(k,j)*~index~(V,Z+j), j= 0..k))
                )[]},
                k= 0..K
            )
        ),
        ICs
    )    
end proc
:

In order to apply this to your system, I needed to make up 6 initial conditions. You could probably come up with better initial conditions than I could. 

Sys:= {
    x[n]=1/3*(2*x[n-1]*y[n-1]+4*x[n-1]*z[n-1])+1/12*(2*x[n-2]*y[n-2]+4*x[n-2]*z[n-2]),
    y[n]=1/3*(1/4*x[n-1]*z[n-1]+y[n-1])+1/12*(1/4*x[n-2]*z[n-2]+y[n-2]),
    z[n]=1/3*(x[n-1]*z[n-1]+2*y[n-1]*z[n-1])+1/12*(x[n-2]*z[n-2]+2*y[n-2]*z[n-2])
}:
Dsys:= RecToDsys(Sys, [x,y,z], n, {x[0]=-1, x[1]=0, y[0]=1, y[1]=0, z[0]=-1, z[1]=-2});
Dsol:= dsolve(Dsys, numeric, method= classical[foreuler], stepsize= 1):
Dsol(10);
plots:-odeplot(Dsol, [x,y,z](n), n= 0..10);

For the initial conditions that I tried, the solutions either quickly blew up (diverged to infinity) or quickly converged to 0. But you may know the proper initial conditions. My procedure does NOT require that the initial conditions start at index 0.

Simply replace a with a/((max-min)(a)/2). For example,

plot(<t | a/((max-min)(a)/2) >);

Use parentheses, not square brackets. 

Your equations are autonomous, which means that the independent variable n only appears as an index. And your equations are second order, meaning that each iteration can be computed strictly from the previous 2 iterations (and the independent variable if they weren't also autonomous). Just as with systems of differential equations, your system can be converted to an equivalent first-order system of 6 equations. At that point, you have an autonomous system of first-order recurrences, so this Post that I wrote 1-1/2 years ago should help you: "Autonomous first-order nonlinear recurrences with parameters and multiple dependent variables".

While you're studying that Post, I'll show how to convert your system to first order.

@Thomas Dean There should be a help page documenting the distributive properties of all operators over all other operators, and which happen automatically. There are two "partially invisble" operators:

1. Function application: This can be thought of as a binary infix operator (without an apparent symbol) that takes a single left operand and an expression sequence (the contents of the parentheses) as its right operand. In prefix form, this operator is `?()`.
2. Index application: Works the same as function application except that its right operand is the expression sequence in the [...]. Its prefix form is `?[]`.

And there are four "outfix" operators (i.e., the operators surround their operands):

1. list constructor:  [...]: Its prefix form is `[]`. But if the [] have any left operand, it's index application, not list construction!
2. set constructor: {...}: Its prefix form is `{}`.
3. rtable row constructor: <...|...>: Its prefix form is `<|>`.
4. rtable column constructor: <..., ...>: Its prefix form is `<,>`.

And there is a "meta" postfix operator: the elementwise operator ~: Its prefix form is `~`[...], where the [...] contains the prefix form of the operator that is the (left) operand of ~.

All of these are at least somewhat overloadable\overrideable, and that's quite commonly done for the rtable constructors.

I think that index application has the highest precedence of any operator. Function application is very high. If you've read ?evalapply, you now know that function application distributes over almost all other operators.

These help pages are very important for understanding operators:

1. ?operators,precedence (although it doesn't cover the seven operators above)
2. ?use
3. ?evalapply
4. ?object,operators

If you've tried a few solutions with dsolve, you've realized by now that it won't work without a specifying a numeric value of the parameter A. Ideally we'd like to learn how the character of the solution changes with different ranges for A. For this, use the parameters option to dsolve(..., numeric, ...):

restart:
ode:= diff(theta(t),t) = 1+A+(A-1)*cos(theta(t));
Q:= dsolve({ode, theta(0)= -Pi/2}, numeric, parameters= [A]):
PlotQ:= proc(A) Q(parameters= [A]); plots:-odeplot(Q, args[2..]) end proc
:
MultiPlot:= proc(As::list(realcons))
    local k, N:= nops(As)-1; 
    plots:-display(
        seq(PlotQ(As[k+1], t= -2..2, color= COLOR(HUE, .85*k/N)), k= 0..N), 
        _rest
    )
end proc
:
plots:-animate(
    MultiPlot, 
    [['seq'(k, k= K..K+2, .05)], title= 'sprintf'("A = %4.2f..%4.2f", K, K+2)], 
     K= -3..1, frames= 200, paraminfo= false
);