@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.

@samiyare You are right. That was my silly mistake!

@samiyare You can use continuation in N[bt]. You get results when N[bt] is 0.4. Then use continuation from 0.4 to 0.2, i.e. set
N[bt]:=cc*0.2+(1-cc)*0.4; #with cc in 0..1.
then N[bt] starts at cc=0 by being 0.4, where dsolve doesn't have a problem and works its way up to cc=1 where N[bt]=0.2. It happens to work.
You need only insert the optional argument 'continuation'=cc in dsolve.
I upload my version of your latest example.

MaplePrimes13-12-27o.mw

(Note the link should work now)

The error comes from dsolve. In your equation eq1 there is a term diff(u(eta),eta)/eta. You are considering the interval 0..1. Thus there will be a problem at eta = 0.
In dsolve you can specify the method as method=bvp[midrich]  (or middefer).
I didn't have any luck right away though. I got an error "Inintial Newton iteration not converging".
You probably need to do some exploration before embarking on the ambitious project of determining p2 and phi0 with the required properties. So just trying to make dsolve come out with some result for suitable values of p2 and phi0 would be a first step.
Comparing your equations with the equations considered on November 27:
http://www.mapleprimes.com/questions/200329-Automatic-Calculation-In-Dsolve-

I notice that there was no similar singularity. In that case the singularity would have been at eta=1 had it not been for the fact that the boundary condition was given at eta=1-zet, where then zet=0.5

Here you assign to zet twice. First zet:=0, secondly zet:=0.5: That of course just means that zet is the last value, but it does make you wonder. In Q the integrals are over the interval 0..1-zet. So that would be 0..0.5.

@Majmaj The output from print is NULL. That it prints to the screen does not in the usual Maple sense mean that it produces an output. You would have to copy and paste the result from the screen. You cannot save it to a variable.
You may want to refer to ?print where it states:
"The function print displays the values of the expressions appearing as arguments, and returns NULL as the function value. "

Try this:
restart;
generate_y := proc (y)
    y := (rand(1 .. 5))();
    print(y);
end proc;
qwerty:=12345;
generate_y(k);
%;
k;
yy:=generate_y(q);
#Compare with
generate_y2 := proc (y)
    y := (rand(1 .. 5))();
    NULL 
end proc;
zz:=generate_y2(p);

Just a follow-up.
The following extremely simple example shows an important difference between Maple 15 and Maple 16 and 17:

Maple 16 and 17:
B:=RootOf((x-n)*(x-2*n)*(x-3*n)*(x-4*n),x);
asympt(B,n,2);
                      
allvalues(%);
                       

Maple 15:
B:=RootOf((x-n)*(x-2*n)*(x-3*n)*(x-4*n),x);
asympt(B,n,2);
                        n

@John Fredsted At the danger of belaboring the point: The matrix M is only hermitian if the names in the entries are real.
op(M); #Showing shape=[]
type(M,'Matrix'(hermitian)); #Reasonable
LinearAlgebra:-HermitianTranspose(M) - M; #OK too
simplify(%) assuming real;
Only after assuming (as evalc does) that the entries are real do we get the zero matrix.