@ The only Classic version of Maple I have on this computer (64bit) is Maple 8.
I tried a variant of the method I described.
That worked well.
I did:

stopat(s1,4);
s1(2,x0);
Up came the debugger prompt DBG>
I entered
next (followed by enter).
Then
Y:=eval(y);
then
quit
and
eval(Y);
unstopat(s1);

@ This is in
`Standard Worksheet Interface, Maple 2016.2, Windows 10, October 21 2016 Build ID 1174570`

It is a great title. It tells us a whole lot of what the subject is. Besides it is a very polite title.

@Carl Love What I did (in Maple 2015.2 and in Maple 2016.2) was:
restart;
Digits:=15:
kernelopts(floatPi=false); #Inserted explicitly here (is also in my maple.ini)
G := 6.6743*10^(-8); c := 2.99792458*10^10; epsilon := 7.3214*10^35;
sys2 := diff(P(r), r) = -G*(epsilon+P(r))*((8.98*10^14*(1/3))*Pi*r^3+4*Pi*r^3*P(r)/c^2)/(c^2*(r^2-(2*G*r*(8.98*10^14*(1/3))*Pi*(r^3))/c^2)), P(0.1e-9) = 8.5561*10^34;
dsolve({sys2});

I recieved the error described above from Maple 2015.2 and the unstoppable behavior from Maple 2016.2.

Maple 18.02 behaves the same as Maple 2015.2 except for an extra error due to the nonexistence of the floatPi kernel option.

I notice that using RGB works (in Maple 2016.2, but not in Maple 2015.2) as in

color = RGB(sin(x),sin(y),sin(z))

@Carl Love  I have kernelopts(floatPi=false) in my maple.ini and in trying
dsolve({sys2});
I got the following error in Maple 2015.2:

Error, (in signum) too many levels of recursion
Modifying the input with evalf gave an immediate answer:

dsolve(evalf({sys2}));

First I tried in Maple 2016.2 (I don't have 2016.1 anymore) and the command dsolve({sys2}) seemed to never stop evaluating. I had to close Maple. Using evalf as above worked.

In both versions starting with kernelopts(floatPi=true) made the command work without the evalf addition.
I shall report the issue as an SCR about Maple 2016.2.

@umar khan You can use v(0) = v0 (where v0 is unassigned) and proceed as I did.
The line

res2:=dsolve(eval({sys2,sys3},res1),numeric);

  you change to

res2:=dsolve(eval({sys2,sys3, w(688240)=-2.05684},res1),numeric);

Then dsolve will with the extra condition w(688240)=-2.05684  (if possible)  find v0 as well as u and w.

You will need to use method=bvp[midrich] or method=bvp[middefer].

I had no luck with my attemps though.
 

@vv Thanks. That should do it with _MaxSols set sufficiently high.
_MaxSols:=200 should do it.

@vv Indeed, solve does a good job (here for the case m=125):

r:=solve(7/18-(1/2)*cos(15625*Pi*x)=0,x,allsolutions);
S:={op(indets(r,suffixed(_B,integer)))=b,op(indets(r,suffixed(_Z,integer)))=z};
R:=subs(S,r);
Here b = 0 or 1, z an arbitrary integer.
But he only wants roots in the interval 6/125 .. 7/125 (still for m=125).

@  I get into problems before you do.

For the record let me collect the code:
 

restart;
CreaCos := proc (C, n, m, t) local k, F; 
   F := C[1][1]+(C[1][2]-C[1][1])*t; 
   for k to n-1 do 
     F := F, C[k+1][1]+((1/2)*C[k+1][2]-(1/2)*C[k+1][1])*(1-cos(Pi*m^k*t)) 
   end do; 
   return F 
end proc;

CreaC := proc (n, c1, c2) local i, H; 
  H := NULL; 
  for i to n do 
    H := H, [c1, c2] 
  end do; 
  return [H] 
end proc;

Rdensf := proc (C, V0, n, m) local k, k0, R, raiz, V, t1, alpha; 
    alpha := evalf(sqrt(n-1)/m); 
    t1 := -10; 
    V := NULL; 
    for k to m do 
      if evalf((k-1)/m) <= evalf(V0[1]) and evalf(V0[1]) <= evalf(k/m) then 
        k0 := k; 
        break 
      end if 
    end do; 
    R := sort([RootFinding:-Analytic(CreaCos(C, n, m, x)[n]-V0[n], x, re = (k0-1)/m .. k0/m, im = 0 .. .1)]); 
    for raiz in R do 
      if evalf(Im(raiz)) <= 1/100000000 then 
        V := evalf(<CreaCos(C, n, m, raiz)>); 
        if evalf(LinearAlgebra:-Norm(V-V0)) <= alpha then 
          t1 := raiz; 
          break 
        end if 
      end if 
    end do; 
    return t1 
end proc;
##
## I run into problems when m=125:
t0 := Rdensf(CreaC(3, 0, 1), <1/18, 1/9, 1/9>, 3, 50);
t0 := Rdensf(CreaC(3, 0, 1), <1/18, 1/9, 1/9>, 3, 100);
t0 := Rdensf(CreaC(3, 0, 1), <1/18, 1/9, 1/9>, 3, 124);
t0 := Rdensf(CreaC(3, 0, 1), <1/18, 1/9, 1/9>, 3, 125);

I ran the code in Maple 2016, 2015, and 15.
In the latter case my patience ran out for the case m=100, but m=124 was OK.
All three complained at m = 125.

The call to RootFinding:-Analytic for m=125 is
RootFinding:-Analytic(7/18-(1/2)*cos(15625*Pi*x), x, re = 6/125 .. 7/125, im = 0 .. .1);
and that for me produces an error, where the corresponding command for m=124 just computes the roots.

@John Fredsted Apparently MaplePrimes was just updated.
Here is the same test today:
eq := a+b/(s+t)+c = int(f(x), x = a .. b)

So the result of using the context menu in Maple using copy special/copy as image is now that you get the editable code and the image, not quite the intention, but it may have its advantages.

Nobody answered yet. That may be for the simple reason that we cannot see the command that produced the error.
Your sentence
   " Then, for k=50, 100, 150... the instruction "
is not complete. What is the instruction?
Also I notice that you write:  "These are the procedures:"
and we only see one, CreaCos.

@umar khan In your uploaded worksheet you have not made the two changes that I suggested and which Tom Leslie seconded:

The Pi/pi problem and the exp/e problem.

That is OK, but you ought to tell us why not.

@tomleslie student_md did implicitly say that a thing like 2+3 works, so an example of something that doesn't work would be good.

@student_md Correct. Nothing is assigned to x1(t) or x2(t). Although that can be done rather easily, I would advise against it. I find it more convenient to leave x1 and x2 unassigned so they can be reused in another ode (or the same ode) in the same session without unassigning.
Here you can do:

ode2:=diff~(X,t)=A.X+B;
ic:=eval(X,t=0)=<-3,5>;
sol:=dsolve({ode2,ic}); # Giving the result a name
## The following construction may look weird the first time, but it ensures that X1 is assigned to the x1(t) result and X2 to the x2(t) result:
#
X1,X2:=op(subs(sol,[x1(t),x2(t)])); #Using capital X1, X2
# To understand it, let us split it up:
subs(sol,[x1(t),x2(t)]);
## then remove the list brackets to get a sequence of length 2
op(%);
X1,X2:=%; #Here assignment to both X1 and X2.
#Plotting
plot(X1,t=-2..5);
plot([X1,X2],t=-2..5); #Both
plot([X1,X2,t=-2..5]); # Phase plane plot: Notice the position of ].