@Tamour_Zubair rkf45 works with error toler ances and with no fixed stepsize. The algorithm seeks to assure that the tolerances are satisfied. If it is impossible an error occurs. It is much more complicated than rk4.
Your rk4 code doesn't just need updating: It has to be rewritten.
You can search the internet for code. What you find might very well differ from Maple's implementation.
#######################
An illustration using the simple ode from above where we actually have the exact solution for comparison:
 

restart;
Digits:=15;
ode:=diff(y(t),t)=y(t);
dsolve({ode,y(0)=1});
res_rk4:=dsolve({ode,y(0)=1},numeric,method=classical[rk4],stepsize=0.01,output=listprocedure); # stepsize=0.01
Yrk4:=subs(res_rk4,y(t));
res_rkf45:=dsolve({ode,y(0)=1},numeric,method=rkf45,output=listprocedure, abserr=1e-15, relerr=1e-12);  
Yrkf45:=subs(res_rkf45,y(t));
T:=Vector(11,i->i*0.1,datatype=float);
Valsrk4:=Yrk4~(T);
Valsrkf45:=Yrkf45~(T);
ValsExact:=exp~(T);
Error_rk4:=Valsrk4-ValsExact;
Error_rkf45:=Valsrkf45-ValsExact;

max(abs(Error_rk4)); #2.73097100489395*10^(-10)
max(abs(Error_rkf45)); #5.52891066263328*10^(-13)

@Tamour_Zubair Your code appears quite correct. Maple uses a version where the k's are defined directly as the f-values, but that the result is the same, as is seen here, where only symbols are used:

restart;
## First your version
rk4_k := proc(f, t, y, h) local k1, k2, k3, k4, ynew; 
  k1 := f(t, y)*h; 
  k2 := f(t + 1/2*h, y + 1/2*k1)*h; 
  k3 := f(t + 1/2*h, y + 1/2*k2)*h; 
  k4 := f(t + h, y + k3)*h; 
  ynew := y + 1/6*k1 + 1/3*k2 + 1/3*k3 + 1/6*k4; 
  return ynew; 
end proc;
## New version with Maple's k's (here called K):
rk4_K := proc(f, t, y, h) local K1,K2,K3,K4, ynew; 
K1 := f(t, y); 
K2 := f(t + 1/2*h, y + 1/2*K1*h); 
K3 := f(t + 1/2*h, y + 1/2*K2*h); 
K4 := f(t + h, y + K3*h); 
ynew := y + 1/6*h*(K1 + 2*K2 + 2*K3 + K4); 
  return ynew; 
end proc;
rk4_k(f,t,y,h);
rk4_K(f,t,y,h);
simplify(%%-%); # 0

A simple example will illustrate the algorithm:
 

Digits:=15:
ode:=diff(y(t),t)=y(t);
dsolve({ode,y(0)=1}); # The exponential function as it should be.
# The numerical result from Maple
resM:=dsolve({ode,y(0)=1},numeric,method=classical[rk4],stepsize=0.1);
resM(0);
resM(0.1);
### The procedure f in this case (NOTE: f  in the algorithm is a function of 2 variables t and y).
f:=(t,y)->y;
h:=0.1;
rk4_k(f,0,1,h);
rk4_K(f,0,1,h);
resM(h);
### Now go 10 steps and write the results to a Vector:
Y:=Vector(11,datatype=float); # Maple uses that datatype in dsolve/numeric
Y[1]:=1; # The initial y0 = 1 (t=0)
t:=0:
for n from 1 to 10 do
  Y[n+1]:=rk4_K(f,t,Y[n],h);
  t:=t+h
end do;
Y[11]; # 2.71827974413517
resM(1); # Perfect agreement

@Joseph Poveromo Yes, that makes quite a difference.
Notice that you cannot use square brackets instead of parentheses.

restart;
ODE:=(2/3)*int(diff(y(x),x)*x^2/(x^2 -1),x)^(-3/2) =int(-sqrt(2*y(x)),x);
infolevel[dsolve]:=3:
DEtools:-odeadvisor(ODE);
sol:=dsolve(ODE);

You will find information about the type y=G(x,y') by going to the help page:
?odeadvisor, patterns

@9009134 I suppose that that could be done easily if you are familiar with Excel. I haven't used it for many years.

I do suggest, however, that you export the whole thing in the form of an (N+1)x(N+1) Matrix, where in your case N=49.
The element (1,1) will be of no importance, so give it some value like -1.
The rest of the first row will be X and the rest of the first column will be Y.
The rest of the elements (i,j) (i,j = 2..N+1) will be given by data[3][i-1,j-1].
Here is a way to do that:
 

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
N:=49;
M:=Matrix(N+1,(i,j)->piecewise(i=1 and j=1,-1,i=1 and j>1,X[j-1],j=1 and i>1,Y[i-1],data[3][i-1,j-1]));
### Checking:
M[1,1];
M[8,10];
M[N+1,N+1];
data[3][N,N];
M[1,2..N+1];
X;
M[2..N+1,1];
Y;

Maple has ways to export a Matrix to a file. There is even an Excel add-n. See ?Excel, and or ?Export.

You didn't attach any file! Uploading a file is done from the MaplePrimes editor: Use the fat green up-arrow.

@9009134 To see all elements in data[3] use
interface( rtablesize = 49);  # see under ? interface in the help.
# Try it, but it takes up a lot of space.
Notice that lists like X and Y are indexed from 1 (not 0) and so are matrices.
data[3] is an Array that has index ranges 1..49, 1..49.
The command:
M:=convert(data[3],Matrix);
makes M a Matrix. data[3] remains an Array.

In Maple 2023.1 and in Maple 2022.2 the result is
Int(-1/((1 + LambertW(1 - u))*u), u).

In Maple 2021.2 the result given is the same wrong answer as the one you report.

@KIRAN SAJJAN The procedure used is `dsolve/numeric/bvp`. It doesn't use any finite element method or shooting method. It uses finite difference.
Here is a simple bvp, you can use to see the way the code is proceeding: Uncomment the debug line.
The code itself can be seen using the showstat command.

restart;
odesys:={diff(y(x),x,x)+sin(y(x))=cos(x),y(0)=1,D(y)(1)=0};
#showstat(`dsolve/numeric/bvp`); #To see the long code
#debug(`dsolve/numeric/bvp`);    # To debug
res:=dsolve(odesys,numeric);
plots:-odeplot(res);

The explanation could be that CubeRoot is seen as superfluous, since Mma has a Surd just like Maple's surd.

But it seems that Surd has been overlooked too!

with(MmaTranslator);
FromMma(`Surd[-32,5]`);

@sand15 In Maple 2023 we get:
 

int(0, x=0..+infinity) = limit(int(0, x=0..a), a=+infinity);
## Result 0 = 0

Which is perfectly fine.
Incidentally, when you mention "the definition" then you must be talking about the improper Riemann integral.
There are several other integrals (e.g the Lebesgue integral), but as Carl Love mentions all return 0 on the input 0 on any domain on which they are defined.

A note about using assuming:
 

int(a*x,x=0..2) assuming a=7; # 14
int(C, x=0..+infinity) = limit(int(C, x=0..a), a=+infinity) assuming C=0; #  0=0
int(C, x=0..+infinity) = limit(int(C, x=0..a), a=+infinity) assuming C>0; #  infinity=infinity
int(C, x=0..+infinity) = limit(int(C, x=0..a), a=+infinity) assuming C<0; # -infinity=-infinity

The equality case works as commented in Maple 2021, Maple 2022, and Maple2023 (at least).

@mmcdara To me it is obvious mathematically that any definite integral of 0 should return 0 pure and simple.
It is another question what happens when you first integrate a (symbolic) constant like C over an infinite range.
To answer your question:
 

z :=0:
i := +infinity:
z*i ;  # = ???

Maple 2023 returns undefined which is reasonable to me and I don't consider that inconsistent.
The same answer comes from direct input:
 

0*infinity;

Answer: undefined. Also fine!

In Maple 2023 the two integrals return 0 whether in 1D or 2D.
Not so in Maple 2021 and Maple 2022. I don't have Maple 2019 nor 2020 on this machine.
The help page does mention the syntax for the ranges as required to be in a list. It seems not to make a difference in any case.
 

@Dkunb You are using the the same initial conditions as in the first round. So in order to pass over to the interval 1.1..2.5 the integration must pass the obvious singularity at R=1.
So where do you want to start on the other side of 1 and in case of T with which value?

@Brian Hanley Try your code without using the ColorTools package, but using [gray,gray] instead.
 

restart;

with(plots):
#with(ColorTools);
#cGr4s := Color([0.50, 0.50, 0.50]);

contourplot3d(-5*d/(d^2 + y^2 + 1), d = -3 .. 3, y = -3 .. 3, color = black, thickness = 3, coloring = [gray, gray], contours = [-2, -1.5, -1, -0.5, 0, 0.5, 1, 1.5, 2]);

contourplot3d(-5*d/(d^2 + y^2 + 1), d = -3 .. 3, y = -3 .. 3, filledregions = false, color = black, thickness = 3, coloring = [gray, gray], contours = [-2, -1.5, -1, -0.5, 0, 0.5, 1, 1.5, 2]);

contourplot3d(-5*d/(d^2 + y^2 + 1), d = -3 .. 3, y = -3 .. 3, filledregions = true, color = black, thickness = 3, coloring = [gray,gray], contours = [-2, -1.5, -1, -0.5, 0, 0.5, 1, 1.5, 2]);

For the first two of the contourplot3d commands I get the error:

Error, (in plot3d) unexpected option: coloring = [gray, gray]

The last command works fine.
 
Thus I can't see that the ColorTools package has anything to do with this behavior.

You are assuming that d is nonnegative:
Executing the following doesn't give degrees, minutes, and seconds:
 

dms(-45.2365);

Maybe this would help, it seems so:
 

restart;
dms := proc(d, m := 0, s := 0)
 local con, d1, d2, d3; 
if 0 <> frac(d) then 
   d1 := floor(d); 
   d2 := floor(60*frac(d)); 
   d3 := 60*frac(60*frac(d)); 
   con := cat(d1, `° `, d2, `' `, d3, `"`); 
else 
con := d + m/60 + s/3600; 
end if; 
cat(con,`° `); 
end proc;