Use a table to store the solutions. Here's an example. I have a system of five randomly generated equations in variables v, w, x, y, z with a parameter a. I use a table indexed by a to store the fsolve solutions. Then I plot how x varies with a. (This particular plot is not very meaningful because I did nothing to control which of the system's many solutions was chosen by fsolve.)

eq||(1..5):= 'randpoly([v,w,x,y,z], degree=2)' $ 5:
for a from 0 by .01 to 1 do
     Sols[a]:= fsolve({eq1+a, eq||(2..5)})
end do:
plot([seq([a, eval(x, Sols[a])], a= 0..1, .01)]);

If you need more details, then please post your code.

The last command in the file is

slope[i][j]:= slope(pa[i], pa[j]);

In this command, you are using slope in two different ways---both as a procedure and as a place to store that procedure's results. This causes an infinite loop that eventually causes the kernel to crash.  You need to use a different name on the left side. You could make it SLOPE. Additionally, you may ultimately find it more convenient to use SLOPE[i,j] rather than SLOPE[i][j]. But I can't tell for sure without seeing what you intend to do with SLOPE.

Use RandomTools:-Generate and an `if` expression.

Tosser:= RandomTools:-Generate(float(range= 0..1, method= uniform, digits= 5), makeproc):
Toss:= ()-> `if`(Tosser() < .495, 0, 1):

Two examples of its use:
seq(Toss(), k= 1..10);
                  0, 1, 1, 0, 1, 1, 0, 1, 0, 1
add(Toss(), k= 1..10^5);
                             50325

Yes, the new ?LinearAlgebra,Generic package works over arbitrary commutative rings (many commands in the package require a field). There is no explicit rank command, but it's easy enough to get it from the row echelon form, for which there is a command. There's also SmithForm and BareissAlgorithm.

The help page ?LinearAlgebra,Generic shows an example of matrix computations over GF(3^3).

When dsolve (or any other Maple command) needs to make up a variable name which is not supplied by the user, it begins that name with _. These variables can be used like any other. The only significance of the _ is that it shows that the variable name was added by Maple.

You asked:

Does it mean that the integral can't be evaluated or maybe something else?

There are two separate issues here: the appearance of the _ in the solution and the appearance of an integral in the solution. As to whether the integral can be evaluated: It depends on the color of the integral sign. If it is blue then it can't be evaluated. If it is black, then it may be possible to evaluate it, and you can try to evaluate it using the value command. But, you don't need to remember the color code because you can apply the value command in either case.

DEtools[reduce_order](
     x*(x-1)*diff(y(x),x$2) - x*diff(y(x),x) + y(x) = 0,
     y(x)=x,
     output= solution
);

collect(expand(value(%)), [_C||(1..2)]);

You need to select based on type name rather than type symbol because indexed names such as a[1] are considered names but not symbols.

V:= < seq(`<|>`(selectremove(hastype, eq||k, name)), k= 1..3) >:
V1:= V[..,[1]]:  V2:= V[..,[2]]:
eq:= < < eq1 >, < eq2 >, < eq3 > >:

You need to use LinearAlgebra:-Equal instead of is, because is will not work with Matrices and Vectors.

LinearAlgebra:-Equal(V1+V2, eq);
                              true

You need to take a limit instead of using direct evaluation. Change the last line from

omega[L]:= 9:  t:= 10:  evalf(U[init]);

to the single command

limit(U[init], omega[L]= omega[G]);

and you will get 0.

Another option is to ?overload igcd so that it will take container objects (sets, lists, vectors, etc.) as arguments.

restart:
~set:= x-> convert(x, set):
local igcd:= overload([
     proc(X::~set, $) option overload; :-igcd(X[]) end proc,
     :-igcd
]):

igcd~([ [12,4], {8,2}, <35|5> ]);
                           [4, 2, 5]

See ?overload , ?coercion , ?updates,Maple17,LocalNames , ?dollar,Procedure (section "End-of-parameters marker $").

I used to be annoyed by exactly the same problem. A workaround that I use now is to get rid of prompts entirely. I do this with the command interface(prompt= ""). I have it in my initialization file.

Another benefit of having no prompt is that when I upload a worksheet to MaplePrimes for display (when that is working (*sigh*)), I have an extra two characters per line to work with before MaplePrimes adds a line break (which wreaks havoc with indented code).

@JohnS Picking up where your code left off:

The variables that appear in the solution in the form a = a are the basic (B) variables. The other variables (nonbasic (NB)) are either constant or assigned linear expressions in terms of the basic variables. The first statement separates out the trivial equations of the form a = a.

B,NB:= selectremove(evalb, Sol):

The next statement simply extracts the variables from those trivial equations.

B:= op~(B):
Zs:= [z||(1..3)]:

Now iterate through the basic variables, setting each to 1 and all the others to 0. Call this set of evaluations E. Then evaluate the nonbasic assignments using E. Finally, evaluate the Zs using E plus the evaluated nonbasic assignments.

[seq(
     (E-> eval(Zs, eval(NB, E) union E))({v= 1} union (B minus {v} =~ 0)),
     v in B
)];

So, each sublist above is [z1,z2,z3] evaluated at a particular basic solution.

From a*r^n = b, you get r = (b/a)^(1/n). We can do the problem with symbolic a, b, and n. You can substitute a=3 and b=5 if you want, but n must remain symbolic. If n is not symbolic, then you can't take the limit as n -> infinity. This may be where you went wrong, but it shouldn't have frozen your computer.

restart:
interface(showassumed=0):
assume(a>0, b>0, n::posint):
f:= x-> x^3:
r:= (b/a)^(1/n):

The kth subinterval goes from a*r^(k-1) to a*r^k as k goes from 1 to n.

A:= sum(f(a*r^k)*(a*r^k - a*r^(k-1)), k= 1..n);

limit(A, n= infinity);

Expanding upon nm's Answer, in order for the loop to end, there needs to be a change to p1 or q1 inside the loop. The expression p1-q1 computes something, but doesn't save it; it doesn't change anything. To change p1 or q1, you need an assignment statement p1:= ... or q1:= .... Like this:

GCD:= proc(p::posint, q::posint)
local p1:= p, q1:= q;
     while p1<>q1 do
          if q1 < p1 then  p1:= p1-q1  else  q1:= q1-p1  end if
     end do
end proc:

Change

unapply(g,x);

to

g:= unapply(g,x);

Here is how to make your program much faster. This will also work for values of theta and phi other than 0. I got most plots in about 10 seconds.

The first step is change Digits from 25 to 5. Next, remove the initialization of theta and phi. These will be given values immediately before each plot. In each line with a sum or add, make these three changes:

1. remove the evalf,
2. Change the sum or add to Sum (capital S),
3. Change the upper limit of summation from 50 or 100 to infinity.

In the loop that starts for k from 2 to 5 do change k to kk. Change all occurences of k in the loop!