My guess is that people would find this a more interesting problem if they were told what you are trying to do.
You don't say much about that.
Also I wonder why you are using the linalg package at all. That one was superseded by the LinearAlgebra package with the release of Maple 6 quite a while ago.

@Carl Love I think that the question was raised in a class discussing the stability of the scheme.
In the help page (pdsolve,numeric,methods) it says about the Euler method,
"For PDEs that are even order in space, stability requires a restriction of the form
                                 xord
                          k < a h    
 for some a  depending upon the problem."
So the dramatic difference was "intended" so to speak.

What do you mean by a matrix's selfadjoint matrix?
Given the matrix H, what does that selfadjoint matrix A have to satisfy? What is its relation to H supposed to be?
You can write H (uniquely) as a sum of a hermitian matrix and an antihermitian one very easily:

H = (H+H*)/2 + (H-H*)/2

But is it just (H+H*)/2 you are after?

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.