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

@Joe Riel I only got Iterator.hdb and Iterator.mla in Iterator/lib.
The error message from with(Iterator);  is

Error, invalid input: with expects its 1st argument, pname, to be of type {`module`, package}, but received Iterator

In Maple 17, where the help works I get the same error message.

The version I installed (or thought I installed) is 2.3.4.

@Joe Riel I got the most recent version of the Iterator package from the Maple Application Center.
I have used an earlier version before. I'm using Windows Vista. In Maple 18 the libname seems quite right:
"C:\Program Files\Maple 18\lib", "C:\Users\alsholm\maple\toolbox\Iterator\lib", "."
That User directory is the place where Iterator.hbd and Iterator.mla are located.
However, ?Iterator doesn't bring up the help for Iterator, but something from combstruct. Also with(Iterator); just gives an error message.
In Maple 17, however, where libname still includes the same User directory, ?Iterator does bring up the help page for Iterator. But also here with(Iterator); doesn't work.
Any guess as to what I might be doing wrong?

@nm It is unsafe to use sum since in that case Maple will only execute the procedures a and b twice: once for n=0 and once for general n. The exceptional cases will not be caught by int. If you replace sum with add the procedures a and b will be called maxN + 1 times so that int will work on concrete values of n.
To see this you can insert a print statement in the procedure a:
a:=proc(n) print(n);
           int(f(x)*cos(n*freq*x),x=from_..to_) /denomC;
         end proc;

Your equation is missing two multiplication signs, they have been put in by DJJerome below.

@Chaimongkol I entered:
ode:=(4*Wi^2*t^3+4*t+4*sqrt(1-t^2)*De)*u(t)+(-8*sqrt(1-t^2)*De*t-4*L-4*De^2*t^3)*diff(u(t),t,t)=0;
infolevel[dsolve]:=5:
dsolve(ode);
And got a lot of printout. At one point it said:
-> Trying a Liouvillian solution using Kovacic's algorithm
   <- No Liouvillian solutions exists

Then it continued with something else, but didn't finish with anything before my patience ran out.

@baharm31 Saving w-values to a vector or matrix is easily done. Here is the simple example I used above.

restart;
pde:= diff(w(x, t),x,x)-diff(w(x,t),t)=0;
pi:=evalf(Pi): #Instead of just 3.14.
res:=pdsolve(pde,{w(x,0)=sin(x),w(0,t)=0,w(pi,t)=0},numeric);
res:-plot3d(t=0..1);
W:=subs(res:-value(output=listprocedure),w(x,t));
W(pi,0.1); #test
#To see what is going on m and n are fairly small.
#Increasing x- and t-values along rows and columns respectively.
t0:=0: t1:=1: a:=0: b:=pi: m:=10: n:=8:
h:=(b-a)/(n-1): k:=(t1-t0)/(m-1):
A:=Matrix(m,n,(i,j)->W(a+h*(j-1),t0+k*(i-1)));
A[..,3]; #The third column: fixed x
A[4,..]; #The fourth row: fixed t
plots:-matrixplot(A);