What is t?  What is a? With a unassigned (which it is here) you cannot do

a[t] := a[t],a[0]

Technically, you can make the assignment, in a procedure, however, when later you try to reference a[t] it you get an infinite recursion.

What are a and s doing?  The assignment a[s] := a[s],a[0] will almost certainly give an error.  That is, it is similar to

restart;
x := x+1;
Error, recursive assignment

If the elements are real, a more efficient technique is to first sort them, then select the minimum adjacent difference:

smallestRealDiff := proc(L :: list)
local sL;
    sL := sort(L);
    min(sL[2..] - sL[..-2])
end proc:

It is not clear what you want.  What do mean by "assigning an equation"?  If you want to assign each equation to an indexed variable, then why not do so directly in a do loop:

for i to 3 do
  eq[i] := a[i] = b[i];
end do:

Here's two alternatives.

eqs := {diff(y(x),x,x)+y(x)^2*x^2=x^2}:
ics := {y(0)=y0, D(y)(0)=dy0}:
integ := dsolve(eqs union ics, numeric, 'parameters'=[y0,dy0]):

Various phase-space plots

integ('parameters' = [0,0]):
plots:-odeplot(integ, [y(x),D(y)(x)], 0..10, 'numpoints'=1000);

integ('parameters' = [0,1]):
plots:-odeplot(integ, [y(x),D(y)(x)], 0..10, 'numpoints'=1000);

Alternative that uses 'initial' to reinitialize the integrator

ics := {y(0)=0, D(y)(0)=0}:
integ := dsolve(eqs union ics, numeric):

plots:-odeplot(integ, [y(x),D(y)(x)], 0..10, 'numpoints'=1000);

integ('initial'); # this will tell you list order
integ('initial' = [0,0,1]):
plots:-odeplot(integ, [y(x),D(y)(x)], 0..10, 'numpoints'=1000);

You need to define what you mean by the difference of two lists. Giving the expected output for the following may suffice,

listdiff( [1,2,2,3,3], [2,1] );

Possible answers are 
 1. error (lists must be the same size)
 2. [-1,1,2,3,3] (difference of corresponding elements, augmented with zero)
 3. [2,3,3] (multiset difference)
 4. [3,3]   (erasure?)

Consider doing what I did when you first posted this.  That is, wrap your code into a Maple procedure, then use the Maple debugger to step through the code.  That is, do

myproc := proc()
  (* all your code *)
end proc:

stopat(myproc):
myproc();

That will bring up the Maple debugger.  You can execute each line one at a time, using the next, into, step, etc commands (read the ?debugger help page).

The problem is equivalent to partitioning a set of 2*m objects into m-partitions of 2 objects.  That can be easily done recursively:

part2 := proc(L :: list)
local n,k,p;
    n := nops(L);
    if n = 2 then
        return [[L]];
    end if;
    return [seq(seq([ [L[1], L[k]], p[] ]
                     , p in procname(subsop(1=NULL,k=NULL,L)))
                , k=2..n )]
end proc:

 part2([seq(1..4)]);
            [[[1, 2], [3, 4]], [[1, 3], [2, 4]], [[1, 4], [2, 3]]]

The conversion of the partitioned integers to the desired form is straightforward. 

Another possibility, which allows only polynomials, and not, say, a+sin(b), is

select(type, L, polynom(constant, {a,b,c}));

That will allow a + b as well as Pi, which may or may not be desired.

See ?type,polynom.

Hmm.  Where does t enter the expressions?  Using all numeric values my machine executes that nested loop (P=100, Q=20, R=20) in about 0.75 seconds.  Using evalhf (and declaring the arrays appropriately) it does it in 0.2 seconds.  Leaving all non-arrays as symbolic it took 2.2 seconds.  This is using an i5 processor.

Note that you can use ?seq with the optional step argument to specify the rows:

 A := Matrix(10,2, (i,j)-> cat(a,i,j));          
                                   [a11     a12 ]
                                   [a21     a22 ]
                                   [a31     a32 ]
                                   [a41     a42 ]
                             A :=  [a51     a52 ]
                                   [a61     a62 ]
                                   [a71     a72 ]
                                   [a81     a82 ]
                                   [a91     a92 ]
                                   [a101    a102]

A[[seq(1..10,2)], ..];                
                                 [a11    a12]
                                 [a31    a32]
                                 [a51    a52]
                                 [a71    a72]
                                 [a91    a92]

No, there isn't a good way to do so.  This sounds like a bug.  You might be able to work around it by saving the model, exiting, then restarting MapleSim.  Let me know if that works.

Here's one approach, using ?selectremove to partition the solutions into two lists.  Note that I used ?fsolve rather than solve; that seems appropriate considering the equation contains floating-points (i.e. is not exact).

det := eval(Determinant(eigen), lambda=1);
fsol := [fsolve(det, m, 'complex')];
(re,cx) := selectremove(type, fsol, realcons);

The immediate problem occurs during the first loop:

x4a, x4b := fsolve((x-h)^2+(eq4-k)^2-r^2);

The argument to fsolve evaluates to

(x-2)^2+87.04

which has no real solutions, hence the assignment error.  You could specify the complex option to ?fsolve, but do you want complex solutions?

This is discussed in the ?algsubs help page; see the last paragraph of the second bullet.  You can do

 foldr(algsubs, q_x, rho=rho_an, 1/rho=1/rho_an);
                         / /d   \           /d       \\
                       k |-|-- p| alpha + p |-- alpha||
                         \ \dx  /           \dx      //
                       --------------------------------
                                          2
                                   R alpha