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

1. You have 30 data values for each of S, P1, and P2_P2e. What are the corresponding times?
2. You say that "the initial point" is not known. By that I guess you mean the option initialpoint to LSSolve. By initialpoint in this sense is only meant a start guess just like in Newton's method you need a start guess for the root finding to get started. If not provided to LSSolve it will try to find one itself.
But you write that k1=0.000438, k4=0.0385, keq=2.7385 and T1_p1=36.8 are what you hope to get. Thus you have a fine initialpoint right there. 
3. Are the parameters beta and Tr in K:=t->cos((1/180)*beta*Pi)^(t/Tr) known? If so what are they?
4. It is somewhat strange to me that K as defined above should be just multiplied onto the ode solutions S,P1,P2,P2_e. So K doesn't enter into the ode system but is just applied afterwards? Any explanation for that?
 

@AlexSss If you are used to programming, then you might like to change the default input input from 2D Math Input to Maple Notation (aka 1D Math Input). At least you should try it.
If you do then space will not be interpreted as multiplication and you can write
with  (LinearAlgebra):
with as many spaces as desired.

To change go to  Tools/Options/Display/Input Display/  and change to Maple Notation.
Keep the Output display as it is.
While still in Tools/Options/ go to Interface/Default format for new worksheets. Change from Document to Worksheet.
Now click the button Apply Globally.

Changes take affect when you open a new worksheet.
These changes can be undone by the analogous procedure at any time.
I never use anything but Maple Notation and Worksheet mode.

@Harry Garst I should have mentioned that since the eigenvalue 0 has geometric multiplicity 4, the corresponding 4 eigenvectors seen in the last four columns of V might as well have been any other 4 linearly independent vectors spanning the same 4-dimensional space. They will in any case be orthogonal to the eigenvector representing the nonzero eigenvalue, but they need not be (and are in fact as we see in the symbolic case) not mutually orthogonal (and so not orthonormal). They could be made so by the Gram-Schmidt orthogonalization procedure.
## But using Gram-Schmidt in the symbolic case without first evaluating the symbolic results at the value of your parameters will create huge symbolic results even if you just try the last 4 columns
GramSchmidt([seq(V[..,i],i=2..5)],normalized);
But if you evaluate V at parameters then GramSchmidt is no problem:

 

params:={lambda[1,1] = 0.9,lambda[2,1] =0.8,lambda[3,1] =0.7,lambda[4,1] = 0.85,lambda[5, 1] = .7, theta[1,1] = 0.25,theta[2,2] = 0.21,theta[3,3]= 0.20,theta[4,4] = 0.15,theta[5,5] = 0.35};
V1:=eval(V,params);
GS:=GramSchmidt([seq(V1[..,i],i=1..5)],normalized);
V1n:=Matrix(GS);
fnormal(V1n^%T.V1n); #The identity matrix
OM1:=map(eval,Omega,params);
simplify(fnormal(V1n^%T.OM1.V1n));

@student_md You can save yourselt trouble of this sort by using a slightly different notation for (in this case) X[1] and w[1] in the actual computations. But that may be awkward. Delaying evaluation by using '  ' only does just that. Even in this case because the axis labels are kept in the plot structure and if that is evaluated again, then the assignment to X[1] (or w[1]) shows up instead.
Below I show a device that goes the other way: Changing the name of the label in an invisible way.
Replace w[1] by  ` w`[1] . That is backquote, space, w, backquote.

restart;
w[1]:=67: #Anything will do
plot1:=plot(x^(1/2),x=0..1,labels=[c,typeset('w[1]')]);
plot1; #67 becomes label now
plot1:=plot(x^2,x=0..1,labels=[c,typeset(` w`[1])]); #The little trick
plot2:=plot(x^(1/2),x=0..1,color=blue); #  If the labels are the same no need for them here.
plots:-display(plot1,plot2);

Finally I should mention that by using backquotes you can turn "anything" into a name, as in this silly example, where the name is `56`.
`56`:=123;
 

@sideshow To see what is going in that expression, try the same thing with an unassigned function F.
 

T:=add((D@@i)(F)(0)*m^i/i!,i=0..2);

D is the operator of differentiation of functions (procedures). So D(sin) returns cos.
D(sin)(Pi/2) therefore returns 0 since cos(Pi/2)=0.
(D@@2)(F) is the second derivative of F, so (D@@2)(sin) = -sin, and so (D@@2)(sin)(Pi/2) = -1.
For convenience (D@@0)(F) is simply F (or if you will, the 0th derivative of F).
That is used above in the add command where I have i running from 0 to 2.
## An ugly looking alternative (in my view) is this:
 

Td:=add(eval(diff(F(m),[m$i]),m=0)*m^i/i!,i=0..2);

For concrete functions like the function f defined in my answer above, T and Td return the same:
 

f:=m->tanh(sqrt(q)*z/y+x*m/y);
eval(T,F=f);
eval(Td,F=f);

You may have a look at the SignalProcessing package:
?SignalProcessing
I don't know anything about that, but maybe you will find what you need there.