Note that x and y are congruent mod 1000 if and only if they are congruent mod 2^3 and mod 5^3. 
F_n is congruent to 7 mod 8 iff  n is congruent to 10 mod 12.

F_n is congruent to 124 mod 125 iff n is congruent to 249, 251, 252 or 498 mod 500.

So F_n is congruent to 999 mod 1000 iff n is congruent to 1498 mod 1500.

Do you mean something like this?

> with(LinearAlgebra):
  M := RandomMatrix(2,3,generator=rand(0..5));


                               [4    1    4]
                          M := [           ]
                               [5    1    2]

> subs(1=a,M);
                            [4    a    4]
                            [           ]
                            [5    a    2]

After seeing from the plot where the approximate min and max are:

> fsolve({Ferr,diff(Ferr,u2)},{u2=12 .. 15, fout = 0 .. 1});

{fout = .8267864835e-2, u2 = 13.60485822}

(Maple tag fails again: this should be

{fout = .8267864835e-2, u2 = 13.60485822}

> fsolve({Ferr,diff(Ferr,u2)},{u2=0 .. 3, fout = -1 .. 1});

{fout = -.5947169278e-1, u2 = 1.191135588}

{fout = -.5947169278e-1, u2 = 1.191135588}

If Q is an expression,

indets(Q, And(function, Not (typefunc(anything,mathfunc))))
 

will contain all function calls that Maple does not recognize as mathematical
functions.  In particular,

> Q := (x+1)(x-1);
   indets(Q, And(function, Not (typefunc(anything,mathfunc))));

       {x(x-1)}

However, this won't help you with something like 3(x-1), which is just evaluated as 3.

1) Your second assume wiped out the first assume, but I suppose you want x to be a nonnegative integer.  Nevertheless, it's best to ignore that at first, get the solution for real x, and then look at nearby integers.
2) There is no such thing as IF in Maple.  There is if, and there is  `if`,  but in this case I think what you want is piecewise. 
3) You might note that. you will always have a < 0 and b < 0 if x >= 0, and then c and d are constant.
4) What is the objective function?  There is nothing labelled as "wage cost".

I assume you mean something like this:

> rsolve(y(t) = sum(b(j), j=1..y(t-1)), y(t));

It's unlikely that rsolve would succeed on this.  rsolve doesn't handle even very simple nonlinear equations very well.  Actually, I don't think there are closed form solutions except in very special cases, or good algorithms for finding them when
they do exist.

The command to plot a surface from data is surfdata in the plots package.  In your case you'll want to use coords=cylindrical.  Thus, if your data form an m x n grid of points:

Data := readdata(myfile, 3):

(result is a list of m*n [r, theta, z] points)

Mesh := [seq([seq(Data[i*n+j],j=1..n)],i=0..m-1)]:
with(plots):
surfdata(Mesh, coords=cylindrical);

The simplest way to do this, I think, is with coordplot in the plots package.

The problem is that it's not enough to produce the plot structure, you must "print" it.
Since you put a colon after your od, that is not happening in your loop.  A cure is to
write

print(plot(...))

instead of just

plot(...)

I downloaded the package, extracted everything from the .zip file to c:\maps, opened the example worksheet in Maple 10 (Standard), and changed the inputs involving file paths:

> libname := libname, "C:/maps";

> load("C:/maps/coarse.m");

It seems to work fine.

You should be using the coefficients of (x-4)^n in fx*qx - px, not the coefficients of
x^n.  You could use a change of variables t = x-4, or series(fx*qx-px, x=4, ...).

You can plot the direction field for a differential equation using DEplot (or phaseportrait or dfieldplot) in the DEtools package.

Maple is correct, and so are you.  In characteristic 7, 3x+5 = 3x+12 = 3(x+4).
The gcd of two members of a Euclidean ring is only unique up to multiplication by a unit (i.e. an element that has an inverse), and 3 is a unit (because 3*5=1 mod 7).

There are several problems with your procedure.

1) You have both a formal parameter f and a local variable f.
Get rid of the local variable.

2) You need f(x)/D(f)(x) instead of f/D(f).

3) It's probably a good idea to use floating point throughout on this.
I've seen students crash Maple in a Newton-Raphson loop, getting rational numbers
with huge numerators and denominators.  So make it

x := evalf(x - f(x)/D(f)(x));

The Matrix command already has a way to define a matrix using a function.  So:

> mat := proc(a,b,c,N)
    Matrix(N,N, (i,j) -> piecewise(j>i+1,a, i>j+1, c, b))
   end proc;

Or you could think of this as a Toeplitz matrix.

> mat := (a,b,c,N) -> LinearAlgebra:-ToeplitzMatrix
       ([a$(N-2), b$3, c$(N-2)]);