The expression  is singular at (x,y)=(0,0):
eval(R0,{x=0,y=0}); ## error
Try another point e.g. (x,y) = (1,0):
 

mtaylor(R0, [x=1, y=0], 5);

I think I managed to read your ode.
I would strongly advocate using odeplot instead of phaseportrait.
You cannot get a direction field (arrows) from odeplot, but your ode is non-autonomous, so a direction field has no meaning anyway.
 

restart;
ode:=diff(y(x),x,x)+8*diff(y(x),x)+9*y(x)+36*y(x)^2=2*sin(x);
## A plot of y(x) versus x:
DEtools:-phaseportrait(ode,y(x),x=-0.2..5,[[y(0)=1,D(y)(0)=0]],numpoints=1000);
# Rewriting ode as a first order system:
ode1:=diff(y(x),x)=v(x);
ode2:=subs(ode1,ode);
## A plot of diff(y(x),x) (i.e. v(x)) versus y(x):
DEtools:-phaseportrait([ode1,ode2],[y(x),v(x)],x=-.2..5,[[y(0)=1,v(0)=0]],numpoints=1000);
#### Much simpler to use dsolve/numeric directly on ode:
res:=dsolve({ode,y(0)=1,D(y)(0)=0},numeric);

## Now plot with odeplot:
plots:-odeplot(res,[y(x),diff(y(x),x)],-.2..5,numpoints=1000,thickness=2);

You may have a look at the help page ?evalf,details.
Also to get inspired you could try showstat(`evalf/constant/Catalan`).
Here is a rather silly example:

restart;
`evalf/constant/K`:=proc()  evalf(sqrt(2)*sin(sqrt(17))) end proc:
evalf(K);
evalf[15](K);
evalf[20](K);
evalf(exp(K)*Pi);
#################
showstat(`evalf/constant/Catalan`);

Since I believe that Maple's solve simply uses a formula for solving polynomials of degree less than 5 the answer that comes out from solve will just be the result of inputting the coefficients of the polynomial into the formula.
Using the very same polynomial as mmcdara, here is a graphical illustration, where the returned sequence from solve is kept in a list, not a set. Thereby we are not relying on any special knowledge of the ordering of sets in Maple.
In the grapical illustration you see that none of the 3 solutions remain real: The first in the list is red, the second blue, and the third is green.

restart;
pol:=x^3+a*x+1;
sol:=[solve](pol,x); # result: a list
with(plots):
complexplot(sol,a=-20..20,color=[red,blue,"DarkGreen"],thickness=3,style=point);
animate(complexplot,[sol,a=-20..A,color=[red,blue,"DarkGreen"],thickness=3,style=point], A=-20..20);

The animation follows   :

You seem to be assuming that (a^b)^c = a^(b*c).
But that isn't correct always: Take a = -1, b=2, and c=1/3:

((-1)^2)^(1/3); 
(-1)^(2/3); 
evalc(%);

a^b ( a <> 0) in complex analysis can be defined as exp(b*log(a)), where log is one of the branches of the logarithm. In Maple log is always chosen as the principal branch (and usually denoted by ln):

exp(2/3*ln(-1));

This isn't a procedure but it is simple:
 

d:=copy(c); # if you don't want to change c
d[..,1]:=floor~(d[..,1]);
d;
c;

If you want to change c:
 

c[..,1]:=floor~(c[..,1]);
c;

You get the animation using odeplot by adding frames=80:
 

p1:=odeplot(Solusix,[x(t),y(t)],0..20,numpoints=1000,color=green,frames=80):
p2:=odeplot(Solusi2,[x(t),y(t)],0..20,numpoints=1000,color=red,frames=80):
p3:=odeplot(Solusi3,[x(t),y(t)],0..20,numpoints=1000,color=blue,frames=80):
display(p1,p2,p3);

You have an error in your ode for n(t):
You have:
 

diff(n(t), t) = alpha_n*(1 - n(t)) - beta_n(t);

It should be:

diff(n(t), t) = alpha_n*(1 - n(t)) - beta_n;

After correcting for that just do:
 

dsolve(dsys1,numeric,method=classical[rk4]);

You can specify the stepsize as in:
 

res:=dsolve(dsys1,numeric,method=classical[rk4],stepsize=0.01);
res(2); #The result at t = 2

Correcting some syntactical problems I got this ODE:
 

ODE:= (2/3)*int(diff(y(x),x)*x^2/(x^2 -1),x) =int(-sqrt(2*y(x)),x);

Is this the one you intended?
If so then the answer in Maple 2023.1 (and also in Maple 2021.2 and Maple 12)  using dsolve is
 

sqrt(y(x)) + 3/4*sqrt(2)/x + 3/4*sqrt(2)*x - c__1 = 0

Differentiating ODE you get a separable ode:
 

ode2:=diff(ODE,x);

From which the resulting implicit solution follows by separation of variables and then integrating both sides wrt. x.

But my interpretation of your original ODE may be wrong?
 

To get the floats do:
 

indets(X1,float);

Save the plot itself in some variable. Below I use p.

p:=%; # saving the plot right after executing the plot3d command.
data:=plottools:-getdata(p);
data[1];
data[2];
data[3];
dx:=5./48; #spacing between x-values
dy:=evalf(b)/48; ##spacing between y-values
X:=[seq(dx*i,i=0..48)]; # List of x-values
Y:=[seq(dy*i,i=0..48)]; # List of y-values
[X[7],Y[9]]; # An example xy point:
data[3][7,9]; # The corresponding value of phi[2]
eval(phi[2],{x=X[7],y=Y[9]}); # Check

As you notice this series is telescoping. So this works:
 

restart;
A:=sum(1/f(n)-1/f(n+1),n=1..infinity); # Using the unassigned f instead of sqrt
B:=subs(infinity=N,A);
expand(B);
eval(%,f=sqrt);
limit(%,N=infinity);

I shall use Maple 12 (since I don't have your Maple 13; later versions have the almost same example in the help for Statistics:-Fit).
B is the independent variable, Bdata is a Vector of values for B, Hdata is a Vector  of values for the dependent variable.

restart;
with(Statistics):
Bdata := Vector([1, 2, 3, 4, 5, 6], datatype = float);
Hdata := Vector([2, 3, 4.8, 10.2, 15.6, 30.9], datatype = float);
H:=Fit(a*B+b*sinh(c*B), Bdata, Hdata, B);

Result: H := 1.28528203195904544*B+.244370171891131000*sinh(.873796260764654664*B)

Plotting H versus B and plotting the data points:
 

plot(H,B=1..6,labels=[B,'H']); p1:=%:
plot(Bdata,Hdata,style=point,symbol=solidcircle,symbolsize=20,color=blue); p2:=%:
plots:-display(p1,p2);

Choosing as initial R, not 0, but here 0.1:
 

restart;
#with(plots):with(plottools):with(DETools):

Sys:=diff(T(R),R)=((1-1/R)*(sqrt(1-(alpha/R)^2*(1-1/R))))^(-1),diff(Phi(R),R)=(alpha/R)^2*(sqrt(1-(alpha/R)^2*(1-1/R)))^(-1);

inits:=[[T(0)=0.5,Phi(0)=0],[T(0)=0.5,Phi(0)=Pi/4]];
K:=dsolve([Sys,op(op(1,inits))],[Phi(R),T(R)],numeric,parameters=[alpha],output=listprocedure);
eval(Sys,R=0);
## Now changing from 0 to 0.1 in inits:
inits2:=[[T(0.1)=0.5,Phi(0.1)=0],[T(0.1)=0.5,Phi(.1)=Pi/4]];
K:=dsolve([Sys,op(op(1,inits2))],[Phi(R),T(R)],numeric,parameters=[alpha],output=listprocedure);
K(parameters=[1]); # Picking alpha = 1 (arbitrarily)
plots:-odeplot(K,[[R,T(R)],[R,Phi(R)]],0..1);

Notice that simplification helps some. It shows that only Phi has a problem at 0, but both have a singularity at 1:
 

simplify([Sys]) assuming R>0, alpha::real;

Output:

Proceeding with Sys2 and changing 0 to 1e-8 in the initialconditions:
 

Digits:=15:
inits3:=[[T(1e-8)=0.5,Phi(1e-8)=0],[T(1e-8)=0.5,Phi(1e-8)=Pi/4]];
K3:=dsolve([op(Sys2),op(op(1,inits3))],[Phi(R),T(R)],numeric,parameters=[alpha],output=listprocedure,abserr=1e-15,relerr=1e-10)
K3(parameters=[1]);
plots:-odeplot(K3,[[R,T(R)],[R,Phi(R)]],0..1);

Here is a version, where A is created as a table:
 

restart;
L:=[seq(2*i - 1, i = 1 .. 10)];
for i in L do
  A[i]:=(x/a)^(i + 1)*(1 - x/a)^2
end do;
eval(A);