Executing your code after correcting a misspelling of MultiSeries (capital S) I got:

asympt(ii_inf,x,3);

Error, (in asympt) unable to compute series
MultiSeries:-asympt(ii_inf,x,1);
Error, (in MultiSeries:-multiseries) unable to sort exponents, {s, 2-s}

Nobody (I'm sure) would want to write the code seen in the image in your link.
So upload a worksheet using the fat green arrow in the MaplePrimes editor.

@Markiyan Hirnyk For your new example I get the exact same output for the two versions:

restart; 
int(ln(-a^2*x^2+1)/sqrt(-x^2+1), x = 0 .. 1);
restart;  
A := Int(ln(-a^2*x^2+1)/sqrt(-x^2+1), x = 0 .. 1):
value(A);

I think this is a limitation and not a bug.
The limitation is also shown in a separated version, where w = u-v and z = u+v. The two resulting pdes for w and z can be solved independently. The pde for w gives no problem, but the one for z shows the same problem as the given system.
 

restart;
sys:={diff(u(x, t), t)-diff(v(x, t), x)+u(x, t)+v(x, t) = (1+t)*x+(x-1)*t^2, diff(v(x, t), t)-diff(u(x, t), x)+u(x, t)+v(x, t) = (1+t)*x*t+(2*x-1)*t};
##
icbs:= {u(0, t) = 0, u(x, 0) = 0, v(0, t) = 0, v(x, 0) = 0};
## 
pde1:=eval(sys[1]-sys[2],u(x,t)=v(x,t)+w(x,t));
solw:=pdsolve(pde1,{w(x,0)=0,w(0,t)=0},numeric,time=t,timestep=0.01,spacestep=0.01,range=0..1);
solw:-plot3d(x = 0 .. 1, t = 0 .. 1); #No problem
pde2:=eval(sys[1]+sys[2],u(x,t)=-v(x,t)+z(x,t));
solz:=pdsolve(pde2,{z(x,0)=0,z(0,t)=0},numeric,time=t,range=0..1,timestep=0.01,spacestep=0.01);
solz:-plot3d( x = 0 ..1, t = 0 .. 1); #Problem

@Markiyan Hirnyk For an equal handling of int and value(Int(...)) clearly assuming should be applied when the integration kicks in.
Here is a more direct version:

restart;
int(ln(-a^2*x^2+1)/sqrt(-x^2+1), x = 0 .. 1) assuming a>0, a<1;
restart;
value(Int(ln(-a^2*x^2+1)/sqrt(-x^2+1), x = 0 .. 1)) assuming a>0, a<1;

@Markiyan Hirnyk Clearly Kitonum is quite right.
Your example certainly doesn't prove he is wrong.
The last part should have been

A := Int(ln(-a^2*x^2+1)/sqrt(-x^2+1), x = 0 .. 1);
res := value(A) assuming a>0, a<1;
which returns unevaluated as did this first one:
int(ln(-a^2*x^2+1)/sqrt(-x^2+1), x = 0 .. 1) assuming a>0, a<1;

@DavidABurton Yes, I can confirm that on Maple 2015.2 and using Windows 10 the plot looks as you show.
On the same machine I tried Maple 18 and Maple 2016.2. In both cases the plots looked much better than the Maple 2015 version, but the Maple 2016 plot didn't look quite as nice (in some sense) as the Maple 18 version.
I also looked at an earlier version, Maple 15. No problem there either.

Maybe one of the experts in these matters can tell us what is going on in Maple 2015.2 w.r.t. graphics.
I suspect that it is a "well known" problem.
 

Since you have tried several things you should upload the worksheet so somebody can look at it.
Use the fat green arrow in the editor.

@nm Here is an non-numeric version:

restart;
ode:=diff(y(x),x$2)+lambda(x)*y(x)=0,diff(lambda(x),x)=0;
bc:=y(0)=0,y(L)=0,D(y)(0)=y1;
res:=dsolve({ode,bc});

Since _2*Z1+_B1 can take on all integer values you can continue with:

subs(2*_Z1+_B1=n,res);
simplify(%) assuming L>0,n>0;

where I have also used yhat we may as well assume that n > 0.

@shahid Trying your worksheet with the (not unimportant, but not necessary) change of replacing 'array' with 'Array' in the dsolve command, I have three remarks:
(1) The order of the columns is given by A1[1,1]  (or eval(A1[1,1]) if you use array).
     This doesn't square with these lines of yours:
    

for j to 1+100*N do f[j] := A1[2, 1][j, 1] end do;
for j to 1+100*N do theta[j] := A1[2, 1][j, 2] end do;
for j to 1+100*N do vel[j] := A1[2, 1][j, 4] end do;

This makes f contain the values of eta, theta the values of f(eta), and vel the values of f '' (eta). This, I'm sure, you didn't intend. In fact you find the values of f in the second column and the values of theta in the seventh column. If vel means velocity and that here is f '(eta), then that is found in the third column.

2. I don't see any shooting done at all. (I don't think there is any reason to do that anyway, though).
3. The title of your reply was "Is this a way to find infinity". That confuses me since you only consider N=1, so that eta=0..1. I had N=20.

@9009134 You get the error message mentioned by Rouben Rostamian:

pdsolve(PDE1, bcs1 union bcs2, numeric);
Error, (in pdsolve/numeric) unable to handle elliptic PDEs

whether or not PDE1 is elliptic!

@Rouben Rostamian  Actually, pdsolve only checks the boundary conditions.
To test that:
(1) Make the correction to the boundary conditions that you mentioned.
(2) Change the pde to the hyperbolic pde PDE1 := diff(H(x, y), x, x) - (diff(H(x, y), y, y)):
(3) Try the command pdsolve(PDE1, bcs1 union bcs2, numeric) once more.

@shesia I'm basically ignorant in the area of statistics.
But the residuals are available from sol:-Results().
Directly by
sol:-Results("residuals");
The order in which they appear as the defining formula for resid in QM is written may not be good for your purpose.
Inside QM use this order instead of the one given there:
 

[seq(seq(resid[i][j],j=1..nops(T)),i=1..3)]

With nops(T)=30 that will make the first 30 rows residuls for P1, the next 30 residuals for P2_P2e, and the last 30 residuals for S.
Here are plots using listplot:
 

plots:-listplot(R[1..nops(T)]);
plots:-display(Array([seq(plots:-listplot(R[(k-1)*nops(T)+1..k*nops(T)]),k=1..3)]));

The first plot:

@Carl Love The code as provided by the OP worked for me in Maple 2015.2 Windows 10, 64 bit.

@shahid I shall try to answer your questions:

1. I used infinity = 4 because that was the value that gave an answer right away.
In order to get a result for larger replacements for infinity you can provide an approximate solution.
If you don't dsolve will supply one itself and that one is based entirely on the boundary conditions.
In this case it is actually [f(eta) = .500000000000000, g(eta) = -0., theta(eta) = -eta+20].
By looking at the graphs obtained when infinity = 4, I came up with this one:
approxsoln=[f(eta)=3*tanh(eta),g(eta)=-sin(2*Pi/inf*eta),theta(eta)=1.2*exp(-eta)]
where inf is the value you will use instead of infinity.
I had no problem with inf:=20.
The comforting fact with the graphs in the case inf=20 is that the curves all 3 level off making the replacement infinity= a finite value like 20 seem reasonable.
When plotting I would use a smaller interval though (0..10, e.g.).
2. I'm not sure that I understand the question.
3. Use the approximate solution above and continuation in epsilon. Code provided below.
4. The options to plot can be used in odeplot as well.
5. Douglas Meade has a package. I have given a method using the parameters option on MaplePrimes some time ago, see
http://www.mapleprimes.com/questions/213212-How-To-Solve-BVP-By-Shooting-Method
################################################
Now the code:

restart;
Digits:=15;
eqn1 := (R/(R-theta(eta))+Omega)*(diff(f(eta), eta$3))+f(eta)*(diff(f(eta), eta$2))+R*(diff(f(eta), eta$2))*(diff(theta(eta), eta))/(R-theta(eta))^2+Omega*(diff(g(eta), eta))+lambda*theta(eta)*cos(alpha)-M*(diff(f(eta), eta))^2 = 0;
eqn2 := (R/(R-theta(eta))+(1/2)*Omega)*(diff(g(eta), eta$2))-2*Omega*(2*g(eta)+diff(f(eta), eta$2))+(diff(f(eta), eta))*g(eta)+(diff(g(eta), eta))*f(eta)+R*(diff(g(eta), eta))*(diff(theta(eta), eta))/(R-theta(eta))^2 = 0;
eqn3 := (1+epsilon*theta(eta))*(diff(theta(eta), eta$2))+epsilon*(diff(theta(eta), eta))^2+Pr*(f(eta)*(diff(theta(eta), eta))-(diff(f(eta), eta))*theta(eta))+Q*theta(eta)+L*exp(-eta) = 0;
#### Leaving out assignment to epsilon:
Omega := 2.: M := .5: R := 5: lambda := 20: Pr := 1: Q := .5: L := .5: W := .5: n := .1: alpha := (1/6)*Pi:
#### The values of epsilon we shall examine:
E:=[ 0, 0.2, 0.5, 1.0, 1.5, 2.0];
bc := f(0) = W, D(f)(0) = 0, D(f)(inf) = 0, D(theta)(0) = -1, theta(inf) = 0, g(0) = -n*(D@D)(f)(0), g(inf) = 0;
#### Using infinity = inf and inf = 20:
inf:=20;
#### The approximate solution found by looking at the graphs when inf=4 and epsilon =0.2 :
AP:=[f(eta)=3*tanh(eta),g(eta)=-sin(2*Pi/inf*eta),theta(eta)=1.2*exp(-eta)];
####
for i from 1 to nops(E) do
   res[i]:=dsolve(eval({eqn1,eqn2,eqn3,bc},epsilon=c*E[i]),numeric,method=bvp[midrich],
              approxsoln=AP,continuation=c)
end do;
#### Since odeplot and display are usede a lot we do:
with(plots):
#### We make an elaborate version only for f:
display(seq(odeplot(res[i],[eta,f(eta)],0..10,color=blue),i=1..6)); pf:=%:
display(seq(display(pf,odeplot(res[i],[eta,f(eta)],0..10,thickness=3,caption=typeset(epsilon = E[i]))),i=1..6),insequence);
#### You can remove 'insequence' below if you don't want an animation:
display(seq(odeplot(res[i],[eta,g(eta)],0..10),i=1..6),insequence);
display(seq(odeplot(res[i],[eta,theta(eta)],0..10),i=1..6),insequence);

Here is the animation of f with all 6 graphs in the backgrund: