@felixiao I copied the code and removed evalm.

I got this:
MaplePrimes23-09-19_long_loop.mw

@felixiao Could you please upload a worksheet: Use the fat green up-arrow in the editor.

You seem to have given us all the code in your question, but I tried copying that and run it (without any evalm).
There were several errors; maybe do the copying and pasting.
So please upload the worksheet.

@Tamour_Zubair It is understandable to me that you want to know about the basics of the numerical methods used to solve odes. But I would leave the rest to the expert(s), who wrote the code for (in this case) Maple's dsolve/numeric. The implementation is superb.
If you see going further as a challenge then OK. Go ahead.

@Carl Love and @Tamour_Zubair

Using the corrections mentioned by Carl I adjusted the code in the worksheet attached.
I made tolerance an optional keyword parameter. The result is returned as a Matrix whose first and second columns consist  of the t and y values, respectively.
### Since the first version I added another version, which in my view is better.
The old one remains and is at the bottom.
the new version doesn't  use lists and allows for use with odeplot.

MaplePrimes23-10-07_rkf45.mw

@Tamour_Zubair The formal parameters used in a procedure cannot be changed inside a procedure:
Let us take a very simple example:
 

p:=proc(h) h:=h/2 end proc;
p(8);

The formal parameter in the procedure p is h. When you apply p to 8 as in p(8), the first that happens is that h is replaced all over the body of the procedure by the number  8. So you have the nonsens assignment 8:=4.
Thus the error message that fortunately reveals the error.
The solution could be this revised version:

p:=proc(h1) local h;
  h:=h1;
  h:=h/2 
end proc;
p(8); # 4

This you can do in RKF45 where the same happens.
 

RKF45 := proc(f, y0, t0, tn, h1) local k1, k2, k3, k4, k5, k6, y, t, Ynext, h_new, tol, ttt1,h; 
  h:=h1;
##The rest unchanged
end proc;

You now have a new problem, however. The computation appears to take forever.

Who wrote the code?

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