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:

@shesia Some of the coefficients in ode_system are given to begin with. That may be justified for some reason. It could be that they are already well established on other grounds.
I'm talking about the fact that we have

k2 := 10^5; T1_s := 14; T1_p2_e := 35; T1_p2 := T1_p1

Accepting the latter one T1_p2 := T1_p1 I tried including k2, T1_s, T1_p2_e with the other parameters T1_p1, k1, keq, k4.
It is not surprising that better results are found when more parameters are allowed to vary.
Results are still not great.
I found
[0.0316135884282020, 0.0128677728179134, 0.209998937003107, 2.18005529815737, 99981.7854333907, 25.5791512364297, 12.0424873109790]
with the ordering [T1_p1, k1, keq, k4,k2,T1_s,T1_p2_e].
The graph of P2(t)+P2_e(t) is (much) better, but the graph of S is worse.

Which version of Maple are you using?
I just tried in Maple 2016, 2015, 18, 15, 12 and got error messages, yes, but not the one you report.
I'm not knowledgeable enough in these matters to say if this is a reasonable problem to begin with.
Next time use text or upload a worksheet so we don't have to type.
The code, if someone else should care to try:
 

eq2:=diff(f(x,y),y,y)-diff(1/x*diff(x*f(x,y),x),x)=200*b*Dirac(y)*Dirac(0.1-x);
bc2:=f(x,0)=0,D[2](f)(x,0)=-w*r0^2/a*Dirac(x-1),f(0,y)=0,D[2](f)(0,y)=0;
pdsolve({eq2,bc2});

@shesia You just gave me 51 time values, but the data for S,P1,P2_P2e only has 30!