Just a comment: The boundary value problem as stated is solved without problem by dsolve.

@J4James Have you tried?
Take e.g. x=0.1. P2[0.1] is a difference between 2 procedures:
P2[0.1];
returns
proc(y) ... end proc-proc(y) ... end proc
Try plotting
plot(P2[0.1],-1-cos(0.1+1)..1+cos(0.1),view=-2..0);
You will see that P2[0.1] is constant as I said above, but it is a procedure in y.
You need the value of the constant. You could use any y in the interval, but the interval is specific to each x. y=0, however, is in the interval for all x.
When plotting a matrix using plot you need a matrix whose columns have numbers as elements, not procedures.
######################
Comment.
I have assumed that your equation is a simple example of something more complicated. If not you could easily solve the problem symbolically.
After having defined d1 and P1 do

res:=dsolve(d1);
res2:=eval(P1,res);
plot(res2,x=0..1);
Notice that the plot looks quite like the plot from the numeric approach.

@Alejandro Jakubi In ?solve,details it states under the heading 'Description':
"The solve command solves one or more equations or inequalities for the specified unknowns. The unknowns may be names, including indexed names (though for efficiency reasons, indexed names should be avoided when possible), or functions. Both indexed names and functions are considered to be independent of each other and of all other unknowns."

So functions ought to be OK as unknowns.

The following two versions do work, though:
p:=(x-3)/(9*x^5-48*x^4+73*x^3-6*x^2-36*x-8);
plot(p, x = -3 .. 3, y = -2 .. 2, discont = [usefdiscont=true]);
plot(unapply(p,x), -3 .. 3, y = -2 .. 2, discont = true);
and the command
discont(p,x);
works fine.
Also doing:
debug(discont);
plot(p, x = -3 .. 3, y = -2 .. 2, discont=true);
shows no sign of a problem although the plot doesn't conform to the findings of discont.

@Kitonum Having defined G already, you can just do:
H:=G@G;

@Carl Love Actually, the square brackets are unnecessary due to the special evaluation rules for evalf.
Notice the difference between the following two calls to evalf:
{x = 4/3+(1/3)*sqrt(7)}, {x = 4/3-(1/3)*sqrt(7)};
evalf(%);
evalf({x = 4/3+(1/3)*sqrt(7)}, {x = 4/3-(1/3)*sqrt(7)});

@mahmood180 If you want to keep the intervals, then don't simplify. You use fnormal and identify so you end up with a piecewise expression, which exactly simplies to the one in which only 0 and 1 are dividing points (only in the simple case in your second worksheet).
To see that clearly try this:
f:=eval(g1(t), sol);
nops(f);
[seq(`+`(op(i..2+i,f)),i=1..7,3)];#List with 3 elements
simplify~(%); #Simplifying each of the 3.
`+`(op(%)); #From list to sum
simplify(%);

@J4James I corrected an error in my reply on singularity. For f to have a fourth derivative (also) at r=-k you need f'''(-k)=0 (in addition to f'(-k)=1), not the equality of f'''(-k) and f''(-k).
This change doesn't really alter anything. The conclusion is the same.
You could try
dsolve(ode-1,f(r),formal_series); #You need the homogeneous eq. But f(r) is a particular solution to ode.
Notice that the two solutions are second degree polynomials and that they are not really different because of the arbitrary constants.
Such a solution will not solve your entire problem because it only has 2 arbitrary constants, so you won't be able to satisfy your 4 boundary conditions.
Another comment:
You can piece together solutions from the intervals -h..-k and -k..h like this:
ode as before.
bcs:=f(-h) = 1/2, D(f)(-h) = 1, f(h) = -1/2, D(f)(h) = 1:
h0:=evalf(1+0.5*sin(1)):
res:=dsolve(eval({ode,f(h) = -1/2, D(f)(h) = 1,f(-k)=h-k+1/2,D(f)(-k)=1},{k=0.5,h=h0}),numeric,method=bvp[midrich]); #Interval -k..h
plots:-odeplot(res,[seq([r,diff(f(r),[r$i])],i=0..3)],-0.5..h0,color=[red,blue,green,black],thickness=[1,1,2,2]); p1:=%:
#And choosing the straight line solution on -h..-k:
p2:=plot([r+h0+0.5,1,0,0],r=-h0..-0.5,legend=["f","f'","f''","f'''"],color=[red,blue,green,black],thickness=[1,1,2,2]):
plots:-display(p1,p2);
You notice that f and f' are continuous,but that f'' is not.
You may experiment to look for other solutions also satisfying the continuity conditions for f and f'.
I haven't found any other.

@mahmood180 I don't understand what you mean.
You could try
u0:=simplify((g2(t)-g1(t))^2); #Giving you an expression piecewise in t
u:=int(u0,t=0..2); #Giving you an expression in the c's
Now minimizing u gives the minimum of u as 0 and the values of the c's as 1.
That is easily checked independently:
eval(u,indets(u,name)=~1); #Returning zero
Clearly sqrt(u) then also has minimal value 0 at those values for the c's.

Try minimizing u:=int((g2(t)-g1(t))^2, t=0.. 2). Clearly that one is minimal where sqrt(u) is.

@J4James Using the simplification suggested by Thomas Richard:
Eq1:=simplify(Eq1,size);
#and then further
ode:=select(has,Eq1,diff);
##
We see that the value of B is irrelevant, but now the singularity at r=-k becomes more visible.
You are trying to solve a bvp on an interval -h..h which has this singularity as an interior point: In your case k=0.5 and h >= 1 (assuming that you always have your epsilon>=0).
No wonder numerical solution is difficult!
In my answer I had k=h=1, thus the singularity was at the left end point forcing me to use a midpoint method.

###
Notice that your equation doesn't involve f, but only the derivatives of f.
ode2:=subs(diff(f(r),r)=g(r),ode);
Boundary conditions g(-h)=1,g(h)=1. But you have the singularity at r=-k of course.
By doing
collect(ode2/(r+k)^3,diff);
you see that if g is going to be 3 times differentiable on -h..h then you need g(-k)=1.
(Edit begin) Also you need
u:=(k+r-1)*(diff(g(r), r, r))+(diff(g(r), r))+(-g(r)+g(-k))/(r+k); #Notice 1 replaced by g(-k)
to have the limit 0 for r->-k. This simply means that the second derivative of g must be zero at r=-k since
limit(u,r=-k);
returns -(D@@2)(g)(-k).
(Edit end)
The bvp for g on -h..-k has the constant solution g(r)=1. On the interval -k..h you have the same solution g(r)=1. Thus g(r)=1 on -h..h is a solution to ode2 so that f(r)=r+C would be a solution to ode. That solution will not solve both of the conditions f(-h)=-1/2 and f(h)=1/2. 
The conclusion is no solution f to the bvp can be 4 times differentiable at r=-k. So you have to find a way to deal with the singularity at r=-k.

@Aakanksha You get the same error in Maple 18. I really don't know anything about MeijerG, but could the parameter lists be inconsistent?
There is a procedure that checks the indices:

showstat(`MeijerG/check_indices`);

From looking at the link
http://en.wikipedia.org/wiki/Meijer_G-function

it appears that the requirements for the indices are that a[k]-b[j] must not be a positive integer for k =1..n, j=1..m.
In your case n=1 and a[1]=-15/2. m=4. So in your case a[1]-b[j] is not a positive integer for any j=1..m since
the b's are -5/2, 5/2, -15/2, -15/2.
Are there other requirements?

In Maple 18 I tried MeijerG on your parameter lists:
param:= [[-15/2], [-13/2]], [[-5/2, 5/2, -15/2, -15/2], []];
MeijerG(param,10.0); #returned .4656987417+0.*I
MeijerG(param,10); #gave the error you reported
So I'm suspecting a bug in the exact evaluation of MeijerG.
I submitted an SCR.

@Aakanksha I tried the code in Maple 9.5 and Maple 10. Both gave results that were quite different from the results from Maple 11 and Maple 18. So something must have changed between Maple 10 and Maple 11.
The plot in Maple 18 (including the holes which are due to small imaginary parts caused by roundoff):

@Aakanksha Which Maple version are you using. In Maple 18 I tried your code (with or without RealDomain) but with Digits at default (10) and it worked fine. There were a few holes in the graph, but using Axel Vogt's suggestion I also tried to plot Re(c) and that graph turned up with no holes.

On this computer I happen to have access also to the much earlier Maple 11. The code worked there too.

restart;
k:=0.99: #with(RealDomain):
m:=1: #Digits:=2:
x:=(Pi*csc(Pi*(k-m)))/(0.693*GAMMA(k)*GAMMA(m));
m1:=MeijerG([[-m],[1-m]],[[0,-m,-m],[k-m]],(-m*k)/snr);
m2:=MeijerG([[-k],[1-k]],[[0,-k,-k],[m-k]],(-m*k)/snr);
c:=x*((((m*k)/snr)^m)*m1-(((m*k)/snr)^k)*m2);
with(plots):
plot(c,snr=0..10);

@belief111 You should take a look at the output of solve(eq,u) in the case you have in your worksheet allvalue.mw, where K:=1. In fact you do that after trying allvalues. You see that solve already gives you all the values although it is using indexed RootOfs. So by doing allvalues(solve(eq,u)) you are applying allvalues to a sequence. That means that allvalues misunderstands the input. Indeed, the 2nd root will be seen as a RootOf-option (see the help) and the 3rd through the 6th roots will be seen as possible other options to allvalues and as such fails miserably.
You may try (still with K:=1)
allvalues~([solve(eq, u)]); #Elementwise application of allvalues
but that doesn't (and shouldn't) give you anything new, but is syntactically correct.
Try after  K:='K'; to execute
solve(eq,u);
you will notice that the output is a single nonindexed RootOf. That represents all the 6 roots and those you can get by using allvalues.