There's no point in making it a function of n and m, if n and m are the dummy variables for the seq.  So just

bcs:=seq(seq(F[n,m](1) = `if`(n=0 and m=0, 1,  0), n=0..nmax), m=0..nmax),
          seq(seq(F[n,m](-1) = `if`(n=0 and m=0, -1,  0), n=0..nmax), m=0..nmax),
          seq(seq(D(F[n,m])(1) = 0, n=0..nmax), m=0..nmax),
          seq(seq(D(F[n,m])(-1) = 0, n=0..nmax), m=0..nmax);
 

You could also take the derivative of f(x) = arctan(x+1) - arctan(x-1)  - arctan(2/x^2), and note that it simplifies to 0, then evaluate, say, f(1) = arctan(2) - arctan(0) - arctan(2) = 0.

You could also take the derivative of f(x) = arctan(x+1) - arctan(x-1)  - arctan(2/x^2), and note that it simplifies to 0, then evaluate, say, f(1) = arctan(2) - arctan(0) - arctan(2) = 0.

This ordering is described in the help page ?RootOf,indexed

This ordering is described in the help page ?RootOf,indexed

> [seq]( v $ min(numboccur(v,a), numboccur(v,b)), v = {op(a), op(b)});

> [seq]( v $ min(numboccur(v,a), numboccur(v,b)), v = {op(a), op(b)});

Though it's not very relevant to Mario's question, the Math 301 prof in me can't resist analyzing Alec's striking picture.  The branch cuts for arctan(z) are (by convention, which Maple follows) taken along the imaginary axis above i and below -i.  So the branch cuts for arctan(z+1) go vertically above -1+i and below -1-i, while those for arctan(z-1) go vertically above 1+i and 1-i.  The branch cuts for arctan(2/z^2) are the inverse images of the branch cuts for arctan under the map z -> 2/z^2, namely the diagonal lines from (+/-) 1 (+/-) i to 0,  Those cuts divide the complex plane into four regions, on each of which arctan(z+1) - arctan(z-1) - arctan(2/z^2) is constant.  Now it's simply a matter of determining the value of that expression at a point in each region: 0 in the black regions and Pi in the red regions. 

Though it's not very relevant to Mario's question, the Math 301 prof in me can't resist analyzing Alec's striking picture.  The branch cuts for arctan(z) are (by convention, which Maple follows) taken along the imaginary axis above i and below -i.  So the branch cuts for arctan(z+1) go vertically above -1+i and below -1-i, while those for arctan(z-1) go vertically above 1+i and 1-i.  The branch cuts for arctan(2/z^2) are the inverse images of the branch cuts for arctan under the map z -> 2/z^2, namely the diagonal lines from (+/-) 1 (+/-) i to 0,  Those cuts divide the complex plane into four regions, on each of which arctan(z+1) - arctan(z-1) - arctan(2/z^2) is constant.  Now it's simply a matter of determining the value of that expression at a point in each region: 0 in the black regions and Pi in the red regions. 

You can specify which root you mean (for a polynomial) by using an index in the RootOf. So RootOf(x^4 - 5*x - 1, x, index = 3) isapproximately -.1996820302 (as you can tell using evalf).  The next problem is that if can't directly compare symbolic constants such as RootOf's.  You can use it in combination with is.  Thus:

for j from 1 to 4 do 
  if is(RootOf(x^4 - 5*x - 1, x, index = j) <= 0) then print(j)
  end if
end do;

           3

(Here root #1 is real and positive, root #3 is real and negative, while roots #2 and #4 are non-real).

You can specify which root you mean (for a polynomial) by using an index in the RootOf. So RootOf(x^4 - 5*x - 1, x, index = 3) isapproximately -.1996820302 (as you can tell using evalf).  The next problem is that if can't directly compare symbolic constants such as RootOf's.  You can use it in combination with is.  Thus:

for j from 1 to 4 do 
  if is(RootOf(x^4 - 5*x - 1, x, index = j) <= 0) then print(j)
  end if
end do;

           3

(Here root #1 is real and positive, root #3 is real and negative, while roots #2 and #4 are non-real).

One way to do it is by expressing the combination of inequalities using max or min.  For example, to plot the solutions of the following system of inequalities :

> ineq[1]:= (x-1)^2 + y^2 <= 3 ;
   ineq[2]:= (x+1)^2 + y^2 <= 4 ;
   ineq[3]:= abs(x) <= y^2;
   combined:= max(seq((lhs-rhs)(ineq[i]),i=1..3)) <= 0;
   with(plots):
   implicitplot(combined, x=-2..2, y=-2..2, grid=[50,50], gridrefine=3, filledregions=true, axes=box);

One way to do it is by expressing the combination of inequalities using max or min.  For example, to plot the solutions of the following system of inequalities :

> ineq[1]:= (x-1)^2 + y^2 <= 3 ;
   ineq[2]:= (x+1)^2 + y^2 <= 4 ;
   ineq[3]:= abs(x) <= y^2;
   combined:= max(seq((lhs-rhs)(ineq[i]),i=1..3)) <= 0;
   with(plots):
   implicitplot(combined, x=-2..2, y=-2..2, grid=[50,50], gridrefine=3, filledregions=true, axes=box);

Cool example.  What's going on is this:

> int(x*sin(Pi*x)*sin(a*x), x=0..1);

Pi*(-2*a+sin(a)*Pi^2-2*cos(a)*a-sin(a)*a^2)/(Pi^4-2*Pi^2*a^2+a^4)

Note that as a -> Pi both numerator and denominator go to 0: this is an "indeterminate form".  This causes big problems for numerical evaluation, which is worse the closer a is to Pi, because roundoff error causes "catastrophic cancellation": both numerator and denominator are evaluated to something very close to 0 but pretty much all roundoff error, and the result is complete garbage rather than the mathematical limit of 1/4. 

As you know, this constant-coefficient linear system should have solutions involving trigonometric functions
(at least for positive values of the parameters).

> sys := [ip(t)+Cp*(diff(vp(t), t)) = 0, Lm*(diff(im(t), t)) = vs(t), Lp*(diff(ip(t), t))+vs(t) = vp(t), ip(t) = im(t)+Cs*(diff(vs(t), t)), im(0) = Ip1, ip(0) = Ip1, vs(0) = 0, vp(0) = 0];

> Try:= {ip(t)=a1*cos(w*t)+b1*sin(w*t),im(t)=a2*cos(w*t)+b2*sin(w*t),
         vp(t)=a3*cos(w*t)+b3*sin(w*t),vs(t)=a4*cos(w*t)+b4*sin(w*t)};

> map(rhs-lhs,eval(sys[1..4],Try));

> eqs:= map(s -> (coeff(s,sin(w*t)),coeff(s,cos(w*t))),%);

> S:=solve(eqs,{a1,a2,a3,a4,b1,b2,b3,b4,w});

S := {a1 = 0, a2 = 0, a3 = 0, a4 = 0, b1 = 0, b2 = 0, b3 = 0, b4 = 0, w = w}, {a1 = -a2*(-1+Cs*Lm*RootOf(1+Lp*_Z^4*Cp*Cs*Lm+(-Lm*Cs-Lp*Cp-Lm*Cp)*_Z^2)^2), a2 = a2, a3 = -b2*(-1+Cs*Lm*RootOf(1+Lp*_Z^4*Cp*Cs*Lm+(-Lm*Cs-Lp*Cp-Lm*Cp)*_Z^2)^2)/(Cp*RootOf(1+Lp*_Z^4*Cp*Cs*Lm+(-Lm*Cs-Lp*Cp-Lm*Cp)*_Z^2)), a4 = Lm*b2*RootOf(1+Lp*_Z^4*Cp*Cs*Lm+(-Lm*Cs-Lp*Cp-Lm*Cp)*_Z^2), b1 = -b2*(-1+Cs*Lm*RootOf(1+Lp*_Z^4*Cp*Cs*Lm+(-Lm*Cs-Lp*Cp-Lm*Cp)*_Z^2)^2), b2 = b2, b3 = a2*(-1+Cs*Lm*RootOf(1+Lp*_Z^4*Cp*Cs*Lm+(-Lm*Cs-Lp*Cp-Lm*Cp)*_Z^2)^2)/(Cp*RootOf(1+Lp*_Z^4*Cp*Cs*Lm+(-Lm*Cs-Lp*Cp-Lm*Cp)*_Z^2)), b4 = -Lm*a2*RootOf(1+Lp*_Z^4*Cp*Cs*Lm+(-Lm*Cs-Lp*Cp-Lm*Cp)*_Z^2), w = RootOf(1+Lp*_Z^4*Cp*Cs*Lm+(-Lm*Cs-Lp*Cp-Lm*Cp)*_Z^2)}

S[2] is the one of interest.

> W:= op(indets(S[2],RootOf));

W := RootOf(1+Lp*_Z^4*Cp*Cs*Lm+(-Lm*Cs-Lp*Cp-Lm*Cp)*_Z^2)

> ws:=allvalues(W);

ws := sqrt(2)*sqrt(Lm*Cs*Lp*Cp*(Lm*Cs+Lp*Cp+Lm*Cp+sqrt(Lm^2*Cs^2-2*Lm*Cs*Lp*Cp+2*Lm^2*Cs*Cp+Lp^2*Cp^2+2*Lp*Cp^2*Lm+Lm^2*Cp^2)))/(2*Lm*Cs*Lp*Cp), -sqrt(2)*sqrt(Lm*Cs*Lp*Cp*(Lm*Cs+Lp*Cp+Lm*Cp+sqrt(Lm^2*Cs^2-2*Lm*Cs*Lp*Cp+2*Lm^2*Cs*Cp+Lp^2*Cp^2+2*Lp*Cp^2*Lm+Lm^2*Cp^2)))/(2*Lm*Cs*Lp*Cp), sqrt(2)*sqrt(Lm*Cs*Lp*Cp*(Lm*Cs+Lp*Cp+Lm*Cp-sqrt(Lm^2*Cs^2-2*Lm*Cs*Lp*Cp+2*Lm^2*Cs*Cp+Lp^2*Cp^2+2*Lp*Cp^2*Lm+Lm^2*Cp^2)))/(2*Lm*Cs*Lp*Cp), -sqrt(2)*sqrt(Lm*Cs*Lp*Cp*(Lm*Cs+Lp*Cp+Lm*Cp-sqrt(Lm^2*Cs^2-2*Lm*Cs*Lp*Cp+2*Lm^2*Cs*Cp+Lp^2*Cp^2+2*Lp*Cp^2*Lm+Lm^2*Cp^2)))/(2*Lm*Cs*Lp*Cp)

The values to use for w are ws[1] and ws[3], as ws[2] and ws[4] are -ws[1] and -ws[3].

> T1:= eval(Try,S[2]);
   Gensol:={seq(v=subs(T1,a2=a21,b2=b21,W=ws[1],v) + subs(T1,a2=a22,b2=b22,W=ws[3],v), v = {im(t),ip(t),vs(t),vp(t)})};

This is the general solution.  Now plug into the initial values and solve for the coefficients.

> consts:=solve(eval(sys[5..8],eval(Gensol,t=0)),{a21,a22,b21,b22});

And finally, the solution:

> Solution:= eval(Gensol,consts);

Well, nobody said it would be simple...