Hint: multiply both sides by one of the denominators

First, you misunderstood the result of listtorec. 

a(n+3)-a(n+2)-a(n+1)+a(n) doesn't mean a(n) = a(n+3)-a(n+2)-a(n+1)+a(n), it means a(n+3)-a(n+2)-a(n+1)+a(n) = 0.  

Now this could be written as

a(n) = -a(n+3)+a(n+2)+a(n+1)

but that would not be useful in your procedure.  Think about why: say if n=3, you're trying to calculate a(3), this formula would involve a(4), a(5) and a(6), and in order to get those you'd need a(7), a(8) and a(9), etc, etc.  See where the "Too many levels of recursion" comes from?

Instead, write the equation as

a(n+3) = a(n+2) + a(n+1) - a(n)

which, with the change of variables n = m-3, becomes

a(m) = a(m-1) + a(m-2) - a(m-3)

and you have something that would be useful...

> plot(sin(x), x=0..2*Pi, legend="Sine of x", 
    legendstyle=[font=[TIMES,ROMAN,12],location=top]);

You mean sin(Pi*n)/(Pi*n).

If n is an integer sin(Pi*n) = 0.  Now for Maple, 0/n evaluates as 0 (unless n has been assigned the value 0, in which case it produces an error). 

In general, for purposes of division a denominator that is not obviously 0 is assumed to be nonzero.  In most mathematical computations this is OK, and it's certainly convenient to not have to worry about whether the denominator could ever be 0 whenever you do a division.  It does, however, cause problems in special cases where the denominator is actually 0.  Look up "specialization problem"

The procedure listtorec in the gfun package will take a list and (attempt to) return a recursion for a sequence starting with that list.  You may need longer lists than what you have.  For example:

> with(gfun):
   listtorec([0,1,4,5,8,9,12,13],a(n));
 

[{a(2) = 4, a(0) = 0, a(1) = 1, a(n+3)-a(n+2)-a(n+1)+a(n)}, ogf]

To check if this is the right recursion:

> P := rectoproc(%[1], a(n)):
   seq(P(n), n=0..20);

0, 1, 4, 5, 8, 9, 12, 13, 16, 17, 20, 21, 24, 25, 28, 29, 32, 33, 36, 37, 40

> plots[display](G, axesfont=[TIMES,ROMAN,12]);

The option for the font on the labels (i.e. the names such as x and y) is labelfont.

The option for the font on the labels of the tickmarks is axesfont.

I don't think it's possible to change the font for just one axis and not the other.  The best you can do is fake one of them using textplot.

First of all, the reason your posting looks strange is that it's in "Plain Text" input format, rather than "Filtered HTML".  For no apparent reason, the Input format sometimes seems to switch to  "Plain Text"  on its own.  Before hitting the "Post comment" button, you should try "Preview comment", and then change the Input format if necessary.

I can't figure out what the problem is because I don't know what you think is unreasonable about the result in file 2031_24.mws.  It looks fine to me.

The error is an interesting bug.

A simpler example:

> limit(E[F], E=infinity);

Error, (in limit/mrv/limsimpl) too many levels of recursion

The point is that E[F] is not really a separate variable independent of the variable of integration E.  If you use EF instead of E[F], and make the assumptions needed to make the integral converge:

> int(E^(3/2)/(exp((E-EF)/(k*T)+1)), E = 0..infinity) assuming k>0, T>0;

(3/4)*exp(EF/(k*T)-1)*sqrt(Pi)/(1/(k*T))^(5/2)

 But the last integral might not have a closed form:

> int(E^(3/2)/(exp((E-EF)/(k*T))+1), E = 0..infinity) assuming k>0, T>0;

int(E^(3/2)/(exp((E-EF)/(k*T))+1), E=0..infinity)

The length of the whole curve y=(x/2)^(3/2) is infinite.  The actual question must be the length of the part of the curve from one point to another.  For the length from x=A to x=B, assuming 0 < A < B:

> Y := (x/2)^(3/2);
   int(sqrt(1+diff(Y,x)^2), x=A..B);

-1/216*(64+18*A)^(3/2)+1/216*(64+18*B)^(3/2)

You didn't say what type of entries these matrices have.  I'll assume they are rational
numbers (in the case of floats, there may be complications due to roundoff). 

One approach is to compute a "signature" for each row, defined in such a way that multiplication of a row by a nonzero number doesn't change the signature.  Rows with the same signature are very likely (but not 100% certain) to be multiples.  For example:

> p := 245428250925877;

(this happens to be prime)

> signature:= proc(V:: Vector)
local n,i,j,T,r;
n:= Dimension(V);
r:= 287998763;
for i from 1 to n do
  r:= 2*r mod p;
  if V[i] <> 0 then
    T:= r;
    for j from i+1 to n do
      r:= 2*r mod p;
      T:= T + r*V[j]/V[i] mod p;
    end do;
    return T  
  fi
end do:
0
end proc;
 

> FindRowMultiples:= proc(M::Matrix)
     local H,i,j,s,r, res;
  use LinearAlgebra in
   res:= NULL;
   for i from 1 to RowDimension(M) do
    r:= Row(M,i);
    s := signature(r);
    if assigned(H[s]) then
      for j in H[s] do
        if Rank(<r,Row(M,j)>) <= 1 then
          res:= res, [j,i]
        end if
      end do;
      H[s]:= [op(H(s)),i];
     else
      H[s]:= [i]
    end if
   end do;
  res
  end use
end proc;
 

Note: this doesn't check for rows of all 0 (such a row is of course a multiple of every other row).

> with(Groebner):

> sys:=  [x1^2-x2^2-x3+x4,
 (x2+1)^2-x5 ,
 125+x1 - x2 ,
 1+x3 - x4 ,
 x8*x6+x7 - x8^2,
 x9*x6+x7 - x9^2-x3 ,
 y2*x6+x7 - x5-y2^2 ,
 y1*x6+x7+x5 - x5-y1^2-x4 ,
 x4+(x1-x2-1)*(x1+x2+1)];

> Solve(sys);

{[[-7813+125*x2, 7812+125*x1, 125*y1^2*x9-125*y1*x9^2-125*x8*y1^2+125*x8*x9^2+125*x8^2*y1-125*x8^2*x9-15751*y1+15876*x9-125*x8, -1968875*y2+15625*y2*x8^2-15625*y2*x9^2-61042969*x8+15625*x8*x9^2-15625*x8^2*x9+63011844*x9+15625*y2^2*x9-15625*y2^2*x8, 63011844*y1-15625*x8^2*y1-1984500*y2-61027344*x8-15625*y2*y1^2+15625*x8*y1^2+15625*y2*x8^2+15625*y2^2*y1-15625*y2^2*x8, -15876+125*x4, -15751-125*x9^2+125*x9*x6+125*x8^2-125*x8*x6, 125*y1*x6+125*x8^2-125*x8*x6-125*y1^2-15876, 15625*y2*x6+15625*x8^2-15625*x8*x6-15625*y2^2-63011844, x8*x6+x7-x8^2, -15751+125*x3, -63011844+15625*x5], plex(x5,x3,x7,x6,x4,y2,y1,x9,x8,x1,x2), {}]}

If you want to see the values for x1, x2, x3, x4, x5 first:

> Basis(sys, plex(x6,x7,x8,x9,y1,y2,x5,x4,x3,x2,x1));

[7812+125*x1, -7813+125*x2, -15751+125*x3, -15876+125*x4, -63011844+15625*x5, -15625*y2*x9^2+15625*y1*x9^2+15625*y2^2*x9-15625*y1^2*x9+61027344*x9-15625*y2^2*y1+15625*y2*y1^2-61042969*y1+15625*y2, -63011844*y1+15625*x8^2*y1+1984500*y2+61027344*x8+15625*y2*y1^2-15625*x8*y1^2-15625*y2*x8^2-15625*y2^2*y1+15625*y2^2*x8, 1968875*y2-15625*y2*x8^2+15625*y2*x9^2+61042969*x8-15625*x8*x9^2+15625*x8^2*x9-63011844*x9-15625*y2^2*x9+15625*y2^2*x8, 15625*y1*x7-15625*y2^2*y1-63011844*y1-15625*x7*y2+15625*y2*y1^2+1984500*y2, 15625*y2*x9^2+1968875*y2-15625*y2^2*x9-63011844*x9+15625*x9*x7-15625*x7*y2, -63011844*x8-15625*y2^2*x8+15625*y2*x8^2-15625*x7*y2+15625*x7*x8, 15625*y2*x6+15625*x7-15625*y2^2-63011844, 125*y1*x6+125*x7-125*y1^2-15876, -15751-125*x9^2+125*x9*x6+125*x7, x8*x6+x7-x8^2]

> solve(%[1..5]);

{x1 = -7812/125, x2 = 7813/125, x3 = 15751/125, x4 = 15876/125, x5 = 63011844/15625}

The <maple> tag seems to be acting up again: that should be

{x1 = -7812/125, x2 = 7813/125, x3 = 15751/125, x4 = 15876/125, x5 = 63011844/15625}

I think you want your first equation to have a right side of LCD1.  But that's not
important.  The fact is that there is no solution for those variables, because there
is an identity not involving those variables:

LEF1*LEF2*Rp^2+LEF1*LEF2*omega^2*LAB1^2-LEF1^2*LAB1*LAB2*omega^2
-LEF1^2*Rp^2+Rs^2*LAB1^2-LAB1*LAB2*Rs^2 = 0

> c,b,a:= seq(coeff(fn,x,j), j=0..2);

In contrast to plot, plot3d is not adaptive, but uses a fixed grid.  You can change the grid option (the default is [25,25]).  So you could try, say,

plot3d(0.2387865921e-3*Pi*sqrt(2)*(1/(Pi*T))^(3/2)*v^2*exp(-0.1924464758e-2*v^2/T), v = 0 .. 1000, T = 0 .. 500, axes = boxed, grid=[100, 100], style=patchcontour);

As for your second question, try

plot([seq(0.2387865921e-3*Pi*sqrt(2)*(1/(Pi*T))^(3/2)*v^2*exp(-0.1924464758e-2*v^2/T),
  T = [100, 200, 300])], v = 0 .. 1000);