constants:= ({constants} minus {gamma})[];

A procedure is a mutable data structure, just like a table, Array, or module.

Use applyrule instead of algsubs:

applyrule(P^q + R_(n+1)^q = B, %);

a) f:= x-> x^3;  g:= x^3;

b) f(2);  eval(g, x= 2);

c) f(y);  eval(g, x= y);

d) x:= 2:  eval(f);  eval(g);

is(5 in {seq(13 &^ k mod 100, k= 1..1000)});

sqrt(`*`(select(isprime, [$1..99])[]));

If you want a for loop to count down instead of up, then you need to include the phrase by -1:

for i from n-1 by -1 to 1 do ...

If you want your output to show what is happening in the loop, then you need to increase printlevel:

printlevel:= 10:
for i from n-1 by -1 to 1 do ...

In your int command, include the option method= _d01ajc. The answer will be returned nearly instantly.

E_clear:=int(I_real, t= t_rise..t_set, numeric=true, method= _d01ajc);
                                       

See ?evalf,int for details.

Replace both lines with the single command

remove(type, identify(d), float);

Matrix multiplication is associative, so grouping in different ways with parentheses is worthless.

If B is invertible, then it is trivial to prove that your expression equals A. So choose a B that is not square but such that B.B' is invertible. (Make B m x n with m < n.)

What your purpose in doing the exercise? That is, Is the goal to perform an exercise of writing the most robust numeric code possible? Or is the goal simply to get the answers? What is the form of the coefficients? Are they exact rational numbers? floating-point numbers? or something more complicated?

If the goal is simply to get the answers, the best option is to use the command solve; it will handle any type of coefficient. Example:

solve(2*x^2 + 3*x + 4 = 0, x);

If you know that you want floating-point answers, then use fsolve:

fsolve(2*x^2 + 3*x + 4 = 0, x);

If your goal is the robust-code-writing exercise, then you need to write a procedure, the skeleton of which looks like this:

SolveQuadratic:= proc(A::float, B::float, C::float)
# Solve A*x^2 + B*x + C = 0 numerically
local ...  #Any other variables in the procedure go here

... #Body of procedure

end proc;

Feel free to ask for more details. All of the above will work in Maple 5.

Here's a start. I have a feeling that this can be improved. I can't test without having your procedure elem.

restart:
Elem:= proc(i,j)
option remember;
     Elem(j,i):= elem(i,j)
end proc:

N:= 2^10:
A:= Matrix(N, N, shape= symmetric);
for i to N do
     A[i,..]:= < Threads:-Seq(Elem(i,j), j= 1..N) >
end do:

Motivated by Axel's analysis, I looked more closely for solutions for the situations that RootFinding:-NextZero rejected. I used a "microscope", searching very close to the origin. By that I mean that I used semilogplot and applied fsolve to numer(f(10^d)) rather than to f(d). I plotted in all cases, and the sequence of plots really shows what's happening.

The first point remains a mystery because the limit(f(d), d=0, right) = -infinity (as Axel noted) where one would expect 0. The limit for the other points is in fact 0, and the plots seem like a continuous progression.

Download erf_solve_2.mw

Download triangle.mw

In order to get a numeric result, you need two numeric limits of integration. Use Int instead of Intat, and then apply evalf.