It looks like you're trying to learn Maple but that you have some experience with other computer languages. The code in my worksheet below may be a bit intense, but perhaps you're already aware of similar things in those other languages. Anyway, please look up what you can in Maple's help pages, and please ask any questions that you have about it. And don't spend too much time trying to understand a confusing help page; you can just ask about it here.

Given that your primary interest in your code is computing the maxes of Collatz sequences, these efficiency improvements are possible:

1. (This one is easy, and obvious): You only need to check for a new max after doing a 3*i+1 step, not after an i/2 step.
2. Memoization (this is subtle): Memoization refers to a procedure recording the results of its calls so that if called again with the same arguments, it doesn't need to redo the computation. Usually, this can be done very simply by putting option remember in the procedure's opening part. That doesn't work very well in this particular case, so my code uses an explicit table as a remember table.

Here's how the remember table works mathematically with the maxes of Collatz sequences (and this would work for any collection of finite sequences of reals sharing a common formula): Suppose that we have two positive integers a and b for which we want to compute the maxes Ma and Mb of their sequences A and B. First we compute A and Ma. We note the position K in the sequence where Ma is attained. In the table, we record that Ma is the max for all sequences beginning with A[k] for k from 1 to K.

Next we compute B and Mb. Suppose that at some point the sequence value in B is one of  those A[k] for which we've already recorded a value. Then it's obvious that Mb is the larger of Ma and the largest value computed so far in B. So, we don't need to finish sequence B! 

Obviously that process (which I just detailed for 2 sequences) can be extended to any number of sequences. If we were only doing a small number of sequences, the extra work required to do this likely wouldn't be worth it. But it is worth it if you want the maxes of a large number of sequences. Below, I do it for all odd numbers upto 1 million in under 7 seconds. Without memoization, it takes over 60 seconds.

Download Collatz.mw

In your 2nd attempt, you are missing a comma after stetsel in the solve command.

In your followup, you didn't load package LinearAlgebra​​​​​. 

The messages should be different; foo doesn't use z. The z-> assign(z) is a separate procedure for which z is a parameter, which is a quite different thing from a local variable. 

I don't know why floor is treated differently than sqrt in this regard---probably an oversight. Anyway, what you want can be done by

plots:-textplot([2,2, Typesetting:-Typeset(floor(5.5))]);

The answer is simply the Taylor series at x=1.

p:= x^4-3*x^2+3:
taylor(p, x=1, degree(p,x)+1);
 

The function that you want is the inverse function of the CDF. In Maple, it's called Statistics:-Quantile. Specifically, you want Quantile(Y, 0.05) for a 5% left tail and Quantile(Y, 0.95) for a 5% right tail.

I do not use the so-called "short-form" names for package commands--ever-- because I think that they reduce code readability (a reader unfamiliar with the packages can't see immediately which package (if any) a command comes from). Instead, at the top level I do something like

LA:= LinearAlgebra:
...
LA:-LinearSolve(...)

Or In a procedure or module:

uses LA= LinearAlgebra;
...
LA:-LinearSolve(...)

I don't think that that example is adequately labelled in the help. I think that it should be labelled "A weird counterexample of what can go wrong due to automatic simplification". If you're just starting to learn Maple, you should ignore this, and probably also the whole use command. I've been a heavy user of Maple daily for over 20 years (with heavy programming experience in several other languages for 20 years before that). I've never used a nested use statement, let alone one whose outer part redefines a variable that is part of the redefinition of the inner part! I say this to assure you that there's no practical reason to learn this example---you'll never need to use it.

The only educational value of this example, in my opinion, is to make you rightfully wary of a somewhat mysterious process known as automatic simplification---which is necessary for efficiency but can lead to unexpected results. It is a preliminary step in the evaluation of expressions that happens before functions are evaluated (the order of evaluation is key to this). The use help page emphasizes that use statements are processed during the automatic simplification. The 2nd part of example 10 should state that the weird result is because b - b is also automatically simplified to 0. If you change that order of evaluation by replacing a - b with its inert form a %- b, and you evaluate that after both use statements are evaluated, then you'll get your expected outcome, (a+b)*a.

If you'd like to know a bit more, I can say more, but I'm only willing to if you have significant experience in Maple or some other language that allows subtle control of the order of evaluation. 
 

The csgn(p) that you got is a result of simplify. It could be avoided (for example, by assume= nonnegative as shown by Ronan), but here's another way. The procedure FractionalExponentAdjuster, below, takes 3 arguments

• P: a pseudo-polynomial, or a set, list, or other structure containing pseudo-polynomials, to be modified;
• x: the main variable of the pseudo-polynomials;
• X: the replacement for x in the returned true polynomials.

It returns two things:

• A modified structure of the same type as P except that the exponents of x have been integerized and x has been replaced with X;
• A procedure that reverts polynomials in X to pseudo-polynomials in x.

FractionalExponentAdjuster:= proc(P, x::name, X::name)
local d:= ilcm((denom@op)~(2, indets(P, identical(x)^rational))[]);
    subs(x= X^d, subsindets(P, identical(x)^rational, p-> X^(d*op(2,p)))),
    P-> subs(X= x^(1/d), subsindets(P, identical(X)^posint, p-> x^(op(2,p)/d)))
end proc
:
#Applying this to your example:
(new,res):= FractionalExponentAdjuster([NSR, deltaa], q, Q);
qu:= quo(new[], Q, 'r');
qr:= res([qu,r]);
#Test. Should get 0 if it's correct:
expand(qr[1]*deltaa + qr[2] - NSR);

Functions that are defined by different algebraic expressions over different real intervals are called piecewise, and the corresponding Maple command is piecewise (see help page ?piecewise). Your function could be entered like this (there are several possible variations to this):

f:= piecewise(Or(x < 1, x > 3), 3*(1-x^2)/((x-1)*(x-3)), 3*x);

Another possibility is 

f:= piecewise(And(x >= 1, x <= 3), 3*x, 3*(x+1)/(3-x));

A crude plot can be made by

plot(f, x= 0..5);

Two ways that it can be improved are limiting the vertical range and telling it to look for discontinuities:

plot(f, x= 0..5, -20..20, discont);

You can see from the plot that f is differentiable at x = 1 because the graph has no sharp point there. Symbolic differentiation of piecewise functions also works:

simplify(diff(f,x));

You can achieve more-accurate results for large x by using 

f:= x-> 4*x*ln(1+1/x);

But the round-off error will still win eventually. So, a better solution is to use an asymptotic series:

f2:= unapply(convert(asympt(f(x), x, 10), polynom), x);

It's fairly easy to customize the number of terms in the series to your value of Digits. The series is alternating, so Leibniz's Theorem on the error of truncated alternating series applies: The absolute error is at most the absolute value of the first omitted term.

Entering small matrices is as simple as this (note that the row separators are semicolons, not commas):

A:= <2, -5, 1; 0, 3, -2>;
B:= <0, 2, -1; 7, 1, 3>;
A+B;

First, remove the epsilon[1] and epsilon[2] from the alias command. Having those be functions of x and t will confuse the next step.  They can be put back as functions afterwards if you need. Then do 

Eq2:= mtaylor(lhs(Eq1), [epsilon[1], epsilon[2]], 2) = 0;

The result can be simplified in various ways. But it's linearized with respect to the epsilons.

You are passing to UP4P a procedure f but UP4P is expecting an algebraic expression f. So, you can correct this immediate problem by changing the last argument in the call to UP4P from f to f(x1,x2,x3,x4).

But, why are you doing this? This program is garbage. You should start from scratch, discard the original program, and try to erase it from your mind (because referring back to this program is corrupting your mind regarding good programming practice).

So, what's your goal?

• learning Maple programming?
• learning general computer programming?
• learning numerical analysis (algorithms)?
• learning numerical analysis (programming)?
• just finding minimum and maximum values of functions?

This program is worthless for any of those purposes.

I don't think that you understand what dualaxisplot is supposed to do. Note that the command is passed two plots (and possibly also some "overall" options such as the title in your case---the extra options are irrelevant to my point here). Look at those individual plots. There's no sense in combining the two plots into a dualaxisplot until the individual plots make sense to you. So, do they?

In both of the individual plots, you're asking it to "plot" three numbers for each value of x. What does that mean? 

Also, x is being treated as a continuous real variable in the interval from 1 to 10, not a discrete integer variable. However, this problem is easy to correct if you can answer the questions from the first two paragraphs.