Actually, it is a bug in the local procedure InputCheck:

showstat(Student:-NumericalAnalysis::InputCheck,1..4);

InputCheck := proc(ODE, IC, trange)
local de, dvar, fde, ivar, range, dvar_i, ivar_r, leftend, inival;
   1   if has(lhs(ODE),'D') then
   2     de := convert(lhs(ODE),'diff') = fde
       else
   3     de := ODE
       end if;
   4   fde := rhs(de);
       ...
end proc

de is defined in line 2 (in your case)  as diff(y(t),t) = fde (not assigned yet).
Then in line 4 fde is assigned the value of rhs(de), which of course is just itself.
Thus fde remains unassigned at this stage and it remains so.
It is as simple as this:

ODE:=D(y)(t) = t*y(t)+1/y(t);
de:=convert(lhs(ODE),diff)=fde;
fde:=rhs(de);
de;

######### I shall report the bug via an SCR.

I interpreted your reply to ThU in http://mapleprimes.com/questions/220473-Getting-Maple-To-Simplify-A-Power-Series   as saying that you would post executable code instead of a picture in the future.
When does the future begin?

If we may assume that gamma is not Euler's constant, but some parameter, then it is pretty simple to find the maximal value for gamma for which the problematic denominator is negative on the whole interval 0..3.
The value is 0.348922854747400.

To get a concise answer you should tell us if you are talking about a function of one or several variables.
Better yet, give us an example of a function of the type you are thinking about.

@umar khan Just copy and paste the code I gave into a fresh worksheet.  You could then just execute, but when copied like that, all commands will be in the same execution group, so all commands will be executed with just one punch on Enter, which is most often not desired.
Use the F3 key on your keyboard to separate into several execution groups, one at a time. Place the cursor in front of a command, press F3 and watch the result. The F4 key does the opposite, it joins execution groups.
I did that on my own computer and the result is this:
MaplePrimes16-12-11odesysCP.mw

This is a minor point, but you use gamma. That is Euler's constant in Maple. Was that intended?
Euler's constant is to 10 digits:
evalf(gamma);
     .5772156649
If I try replacing gamma by 0.1 the problems don't show up.
Using Tom Leslie's notation try

indets(odesys,name);
odesysG:=eval(odesys,gamma=.1);
and then work with odesysG instead.

Thus my question about gamma may not be so minor anyway.
What was gamma supposed to be?

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