@tomleslie Doing

solve(odesys,diff~({F,G,H,thetaP}(eta),eta) union {diff(f(eta), eta$3),diff(theta(eta), eta, eta)});

we find that the denominator in the theta equation is   2.308862660*f(eta)^2-1.
Plotting that denominator using the solution sol1, we see that the denominator becomes zero at a point between eta = 1 and eta = 1.5.
This would certainly give problems:

plots:-odeplot(sol1,[eta,2.308862660*f(eta)^2-1.]);

The corrected syntax is
sol := dsolve(sys union { P(ep) = 1.3668*10^14, m(ep) = 0, rho(ep)=4.2983*10^(14),v(ep)=-0.7342}, numeric,method=rosenbrock, range = 0 .. 20*10^(5));

But you have other problems. I have no time to address those now.

@ My Maple 8 basically only works "under protest". If I try to save anything (or open anything) it closes immediately. I have to run it as administrator each time.
But here is a paste of the screen:

restart;
> s1:=proc(n,x)
> local y,xx,i,j,zz::array(1..n,1..n);
> for i from 1 to n do for j from 1 to n do zz[i,j]:=x[i]*(1+x[j]^2);od:
> od:
> y:=array(1..2,[(1)=x, (2)=zz]):
> for j from 1 to n do xx[i]:=zz[i,i]/(add(zz[i,j],j=1..n));od:
> 0;
> end proc;

  s1 := proc(n, x)
local y, xx, i, j, zz::array(1 .. n, 1 .. n);
    for i to n do
        for j to n do zz[i, j] := x[i]*(1 + x[j]^2) end do
    end do;
    y := array(1 .. 2, [1 = x, 2 = zz]);
    for j to n do xx[i] := zz[i, i]/add(zz[i, j], j = 1 .. n)
    end do;
    0
end proc

> showstat(s1);
>

s1 := proc(n, x)
local y, xx, i, j, zz::array(1 .. n,1 .. n);
   1   for i to n do
   2     for j to n do
   3       zz[i,j] := x[i]*(1+x[j]^2)
         end do
       end do;
   4   y := array(1 .. 2,[1 = x, 2 = zz]);
   5   for j to n do
   6     xx[i] := zz[i,i]/add(zz[i,j],j = 1 .. n)
       end do;
   7   0
end proc

> x0:=Vector(2,[1,1]);
>

                                    [1]
                              x0 := [ ]
                                    [1]

> s1(2,x0);
>

                                  0

> stopat(s1,4);
>

                                 [s1]

> s1(2,x0);
>
2
s1:
   4*  y := array(1 .. 2,[1 = x, 2 = zz]);

> next
array(1 .. 2,[(1)=(Vector[column](2, [1,1], datatype = anything, storage = rectangular, order = Fortran_order, shape = [])),(2)=zz])
s1:
   5   for j to n do
         ...
       end do;

> Y:=eval(y);
s1:
   5   for j to n do
         ...
       end do;

> quit
Warning, computation interrupted

> eval(Y);
>

                              [[1]    ]
                              [[ ], zz]
                              [[1]    ]

> unstopat(s1);
>

                                  []

########## I notice that the debugger prompts DBR> just becomes > when copying and pasting.
Here is a pure text version:
test.txt

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