Taking the max without using abs seems not quite right. So either find max and min or find the max of the absolute value:
errfunc := abs@rhs(errfunc[3]);
Furthermore, since spacestep=0.4 it seems more relevant to use
maxerr := max(seq(errfunc(i*0.4),i=0..5));
and similarly for maxerr2.

I cannot tell if the following idea would be relevant for your situation:
Give the objective function option remember. If computation exceeds the timelimit you still have the computed values available in the memory table. You could sort the entries according to size (keeping track of the indices). You may need to trap error using try ... catch. 

@Markiyan Hirnyk  The limit has nothing to do with n. It is x. I don't think this question deserves more of my time. Sorry!

@Markiyan Hirnyk 

is the result of taking the limit, now a and b are a=(2*j-1)/n*Pi/2, b=a+Pi/n.

It is really not very deep. It uses nothing but Calculus 1 (in the US sense).

tan(n*x) tends to -infinity as x -> a from the right and tends to infinity when x tends to b from the left.
Thus knowing the limits of arctan you have the result.

@Markiyan Hirnyk In the integral the integrand is continuous  and surely limits are taken in the integral resulting in the values for a and b stated, i.e. a=(2*j-1)/n*Pi/2, b=a+Pi/n.

So yes, it is that simple.

@Markiyan Hirnyk Integration by parts:

restart;
A:=Int(u(x)*diff(v(x),x),x=a..b);
                                              
IntegrationTools:-Parts(A,u(x));
                                             

This is certainly valid on any interval [a,b] on which u and v are continuously differentiable. Continuity at a and b is only required of v in the sense that the limits limit(v(x),x=a,right) and limit(v(x),x=b,left) both exist. Similarly for u.
But arctan(2*tan(n*x)) has limit Pi/2 as x->b from the left and limit -Pi/2 as x->a from the right.

@Markiyan Hirnyk The observation that the series is quite slowly convergent for small positive values of s is certainly correct. So the usefulness of the series for those values of s is limited for sure.

@Markiyan Hirnyk Certainly the infinite series diverges at s=0.
S(N) is the Nth partial sum and as such it is finite for all s.

@raveen The exponential function is exp(-I*w*x) not e^(....).

But use fourier from the inttrans package:

restart;
with(inttrans);
k := diff(u(t, x), t,t) = diff(u(t, x), x,x);
bc := u(0, x) = 1/(1+x^2);
v:=D[1](u)(0,x)=0;
eq1:=map(fourier,k,x,w);
eq2:=subs(fourier(u(t, x), x, w)=g(t),eq1);
inits:=g(0)=fourier(rhs(bc),x,w), D(g)(0)=0;
dsolve({eq2,inits});
invfourier(rhs(%),w,x) assuming t::real;

@digerdiga Well, for s <= 0 the integral is divergent.
However, the taylor series about s=1 should have radius of convergence equal to 1.
The coefficients in the taylor series are

(D@@n)(f)(1)/n! = (-1)^(n)*(n+1)*Int(t^2*cosh(2*t)^n/(1+cosh(2*t))^(n+2),t=-infinity..infinity);

The Nth partial sum in the series as a function of N:
S:=N->Sum(A(n)*(s-1)^n,n=0..N);

v10:=value(S(10));
v11:=value(S(11));

plot([v10,v11],s=0..1);

@raveen 
You can do
plot3d(rhs(sol),t=-10..10,x=-10..10);
#More fun with an animation in time:
plots:-animate(plot,[rhs(sol),x=-10..10],t=-10..10);

@Markiyan Hirnyk 
What is happening here:
restart;
with(numapprox);
f:=x->arcsin(sqrt(x));
Digits:=100;
chebyshev(f(x), x=0..1, 0.01);
p:=subs(T=ChebyshevT,%);
evalf(eval(p-f(x),x=0.001)); # Greater than the infnorm result!
infnorm(p-f(x),x=0..1,xmax);
xmax;
evalf(eval(p-f(x),x=xmax)); #Consistent with the infnorm result
plot(abs(f(x)-p),x=0..1);

Note:
The results seem not to get any better with

chebyshev(f(x), x=0..1, 0.0001);

So the question is:
Can uniform approximation with ChebyshevT polynomials for this particular function and on the interval 0..1 get any better than error = 0.016.. ? Maple stops at ChebyshevT(17,2*x-1) no matter what the error is set at in chebyshev.
Experiments with intervals 0.1..0.9, 0.1..1, and 0..0.9 show quite different results. Much higher ChebyshevT polynomials appear.

@Hypnos The RootOf does indeed contain test10. For my example given above I get:


You are not going to be able to do much better since the characteristic polynomial visible inside RootOf is of degree 6:
LinearAlgebra:-CharacteristicPolynomial(A,lambda);
collect(%,lambda);

@raveen 

First about diff versus D.
diff(u(0,x),t) returns 0 because it doesn't depend on t at all.
D(sin) returns cos literally. Thus D takes a function (=procedure) as input and returns a function.
D(sin)(Pi/2)  will return 0 since cos(Pi/2)=0.
The index in the indexed version D[1] refers to differentiation w.r.t. the first variable of a function of several variables. D[1](u) is itself a function of several (in your case two) variables. So D[1](u)(0,x) makes sense.

res:=pdsolve(k); #Getting the general solution
#Now find _F1 and _F2 using bc and v.
eq1:=subs(t=0,bc,res); #with this syntax subs makes its substitutions t=0 and bc in sequence.
#That gave us one equation.
diff(res,t); #Differentiate both sides of res w.r.t. t
convert(%,D); #Want to use v directly so need to convert
eq2:=subs(t=0,v,%); #Again sequential substitution
eq3:=map(int,eq2,x)+(C=0); #We integrate both side thus the need for map (it wasn't needed with diff, but could have been done). 
# ..+ (C=0) added an arbitrary constant C to the left hand side
solve({eq1,eq3},{_F1(x),_F2(x)}); #Solving for the two expressions in the curly brackets
#We now assign _F1 and _F2 to be the functions on the right hand side of the previous result:
_F1,_F2:=op(unapply~(rhs~(%),x));
#Now our res should evaluate to the final result:
res;
#But in order to see that there is no dependence on C and to make the result look nicer we do:
sol:=normal(%);

#Note I omitted an explanation of the assignment to _F1 and _F2. Reason being it is long. But just a few hints: rhs means right hand side. A tilde after any function f (procedure) as in f~ means that f will be applied to each element as in sin~([Pi/2,Pi]) which would correspond to [sin(Pi/2),sin(Pi)].
unapply( x^3+7*sin(x), x) makes the function x->x^3+7*sin(x).
op returns the operands of an expression. For a set or a list the operands are the elements.
Finally a double assignment is made:  a,b:=7,9; means the same as a:=7; b:=9; 

You write that you want to determine the values of y[i]. From what equation? I only see an initial value problem for a differential equation. That ivp is easily solved in closed form as well as numerically using methods Euler or Runge-Kutta or whatever.