The linear system of equations is inconsistent. In particular, the coefficients of tau[1] and tau[4] are identical.

Situations where string manipulations are needed to perform mathematical operations are extremely rare. Your string work can be replaced easily by using eval to change exp, like this:

a:= exp(2*t);
b:= exp(1)^(2*t);
eval(a-b, exp= (x-> exp(1)^x));

Likewise, numeric operations are not needed to verify the equality of expressions; simplify can be used, like this:

answer:= 5*exp(7.5*t);
response:= 5*e^(7.5*t);
simplify(eval(answer-response, e= exp(1)));

or, better yet,

verify(answer, eval(response, e= exp(1)), simplify);

Note that these solutions still work even if the symbol e has not been used in the response. 

The situations that you described as "terrible computations" are unavoidable consequences of the rounding of decimal arithmetic. You'll find them in any system of decimal arithmetic, be it a computer program, a hand calculator, or pencil and paper. It's called loss of significance, or, more dramatically, catastrophic cancellation. You'll find umpteen hits web-searching either of these expressions.

The way to avoid this is to actually use your CAS as something more than a fancy calculator---replace numeric computations with symbolic ones.

A negative number to a fractional (or decimal) power is complex. You have (f '')^0.4.

You didn't try expand(res), which I think will work, possibly followed by simplify(..., symbolic).

Just answering on spec; not at my computer to test it.

Your code is correct, and should not give an error. I suspect an installation or licensing issue.

What do you mean by "multiple roots"? A sequence can only converge to one value. Yes, an iterative root finder can be used to find multiple roots, but only if given multiple initial values, and there'll be separate convergent sequences for each root.

We'll define the absolute error of the kth term of a sequence as the absolute value of the difference between this term and the limit. If a sequence converges, we'll say that the order or convergence is r if the absolute error of the kth term of the sequence can be approximated by A^(r^k) for some 0 < A < 1 and r > 0. The values of A and r can be estimated by linear regression. Here are some Maple procedures for doing that:

restart
:
(*>>>-------------------------------------------------------------------------
This procedure approximates the order of convergence of a convergent sequence
S. It outputs 3 values: the order of convergence, the error-approximating 
function, and the coefficient of determination (aka r-squared) of the 
regression.

This syntax requires Maple 2018 or later.
-------------------------------------------------------------------------<<<*)
ConvergenceOrder:= proc(S::~Vector);
local Err:= (-ln@abs)~(S[..-2] -~ S[-1]), n:= numelems(Err), k, j, R, r;
   Digits:= 15;

   #Check integrity of sequence:
   for k to n do until Err[k] > 0;
   if k > n-2 then
      if k > n then 
         error "sequence doesn't converge sufficiently rapidly"
      fi;
      error "need more terms in sequence"
   fi;
   for j from k to n-1 do 
      if Err[j] >= Err[j+1] then 
         error "absolute errors must be strictly decreasing" 
      fi
   od;

   R:= Statistics:-LinearFit(
      [1, _k], <($k..n)>, ln~(Err), _k, 
      ':-output'= [':-rsquared', ':-parametervector']
   );
   Digits:= 5;
   evalf([(r:= exp(R[2][2])), _e = exp(-exp(R[2][1]))^(r^_k), R[1]])[]
end proc
:
(*>>>----------------------------------------------------------------------
This procedure takes a function f whose root we want and a template T for a
root-finding method that uses only one initial value and returns the 
iteration function.
-----------------------------------------------------------------------<<<*)
MakeIterator:= proc(f::procedure, T::procedure)
local x;
   unapply(simplify(T(x,f)), x)
end proc
:
(*>>>----------------------------------------------------------------------
This procedure applies the iteration J, starting with initial value v0, and
returns the sequence of iterates as an Array.

This uses Maple 2019 syntax.
-----------------------------------------------------------------------<<<*)
Sequence:= proc(
   J::procedure, #unary iteration function 
   x0::complexcons, #initial value
   {  #keyword parameters:
      abserr::And(realcons, positive):= 10.^(1-Digits), #convergence criterion
      max_iters::posint:= 99 #max numbers of iterations
   }
)
local x:= x0, x1, X:= Array([x]), n;
   for n to max_iters while abs((x1:= evalf(J(x))) - x) > abserr do
      X,= (x:= x1)
   od;
   if n > max_iters then error "did not converge" fi;
   X
end proc
: 

The iterative method that you gave is close to---but not quite the same as---Newton's method. If you provide a function and initial value for which your method converges, I'll apply these procedures to it. In lieu of that, I'll just show an example using regular Newton's method:

#Example usage:
Digits:= 999: #Large Digits is needed to get sufficient data.

M:= (x,f)-> (x - f(x)/D(f)(x)): #method--Newton's method in this case.
f:= x-> x^3 - x - 2: #function to find root of
x0:= 1: #initial guess of root

ConvergenceOrder(Sequence(MakeIterator(f,M), x0));
                                /      _k\         
                                \2.0391  /         
            2.0391, _e = 0.77243          , 0.99721

So, the order of convergence is about 2.04, the error-approximating function is as shown, and the coefficient of determination is > 99.7%, indicating that the regression was highly accurate.

I'll give you two methods---the first is easy and lazy and the second guarantees uniform selection from the available primes (meaning each prime is equally likely to be chosen):

EasyRandPrime:= (R::range(realcons))->
   nextprime(rand(ceil(lhs(R))-1 .. prevprime(floor(rhs(R)))-1)())
:
#This syntax requires Maple 2018 or later:
UniformRandPrime:= proc(R::range(realcons))
local rnd:= rand(ceil(lhs(R))..floor(rhs(R))), r;
   do until isprime((r:= rnd()));
   r
end proc:

The usage is the same for both:

EasyRandPrime(3..30);
UniformRandPrime(3..30);

If you require a large number of random selections from the same range, the above procedures can be made more efficient fairly easily.

The Answer by dharr exploits the low-level structure of a procedure in a way that won't work if f's parameters have type decalrations and is potentially dangerous if f's parameters have default values.

The Answer by Christian Wolinski does something that could be done more simply with the commands curry and/or rcurry.

What you want can be achieved much more robustly by simply declaring the procedure f such that its parameters have default values that are global names:

f:= proc(a:= ':-a', b:= ':-b', c:= ':-c') a*b*c end proc:
f(3);
                             3*b*c

Another option, which is even more robust are eliminates the correspondence between parameters and the order that they're declared--but which requires a little more typing--is to use keyword parameters: 

f:= proc({a::algebraic:= ':-a', b::algebraic:= ':-b', c::algebraic:= ':-c'})
   a*b*c 
end proc:
f(b=7);
                     7*a*c

Now f's parameters can be given in any order.
 

The values of keyword parameters of type name and the keywords themselves (the word output in this case) are always global. When there is a conflict between a local (your local Q) and a global (your desired Q setting of the keyword) with obstensibly the same name, the local name takes precedent. In those cases, it's necessary to distinguish the global with the :- prefix, thus, ':-Q'.

I give this Answer just to explain the reasoning behind the other Answer.

Some Maple programmers always use the ':-...' syntax for these unevaluated global names.

It's as simple as 

diff(BesselI(alpha, x), x);

Note that the last letter of the function name is uppercase I ("eye"), not lowercase l ("ell").

I can't get Maple to give the Fourier transform of your function, even by direct integration in the a=1 case.

However, regarding your other question, about algsubs: As discussed on its help page ?algsubs, this command treats denominators differently than non-denominators. To substitute x=ex (where x is a name and ex is an expression) into another expression LE, there are only two commands that should be used: either eval or subs---reserve algsubs for more-complicated substitutions, such as when x is an expression.

• If you want to make the substitution regardless of whether it's mathematically valid, use subs(x= ex, LE).
• If you want Maple to consider the mathematical validity, use eval(LE, x= ex).

Either of these commands will be more robust and far more efficient than algsubs.

For this reason among others, SearchText should be deprecated. StringTools:-Search and StringTools:-SearchAll are externally compiled, so they're adequately efficient replacements; probably more efficient when multiple searches are specified in a single command. Making a wrapper for a case-insensitive version of these to replace searchtext is trivial (if we say that only ASCII characters can have a "case").

While I agree with your premise in theory, in practice you may need to make a distinction between different kinds of infinity. The code that you show returns +infinity, somewhat by happenstance. It could just as well have been -infinity or some complex flavor of infinity. Also note that once you allow infinities, the usual field rules of real and complex arithmetic no longer apply. In general, Maple's ability to navigate these issues on its own is pretty good.

As Acer's Answer shows, mod is complicated, and so I usually reserve its use for "higher" algebraic and number-theoretic computations. For numeric arithmetic, I use irem, which automatically returns unevaluated for non-numeric arguments--suitable for your purpose.

iremp:= (x,m::posint)-> irem(m+irem(x,m),m):
plot(iremp(ceil(x),2), x= -9..9);

The first of these commands is to make the correct adjustment for negative arguments, similar to the distinction between modp and mods.  

A not-necessarily-integer parameter---nominally called "degrees of freedom"---is commonly used for the t and chi-square distributions. See "Welch's t test" and "Welch-Satterthwaite equation". This is taught even in second-semester statistics. Perhaps "degrees of freedom" wasn't the best choice for the name of this parameter, but this usage is unfortunately well established, and the mathematical concept is solid.