@felixiao I couldn't find out what was wrong; everything looked fine.
Therefore I copied the whole code and put it in a text editor (actually Notepad++).
In there I couldn't see any problem either, but I copied the code planted there and put it in a new worksheet in Maple.

Save it as the .mw file:  MaplePrimes23-09-22_Matrix_problem.mw

### Note added later: I did this in Maple2023.1, but I also tried in the much older Maple 12. The problem was the same there. For the remedy this time I used the Notepad that comes with Windows.

####

The plot looks like this:
.

@vv PolynomialSystem and solve both work with V as a list in Maple 12, but not in Maple 2021.

@Carl Love The matrices that are being inverted resulting in H5 and H7 and their determinants:

H51:=Matrix([[G11, G12], [G21, G22]]); ## ~H5
LinearAlgebra:-Determinant(H51);  ## 0. + (9.2517792695382088910002830954*10^(-32))*I
H71:=Matrix([[Z11, Z12], [Z21, Z22]]);   ## ~H7
LinearAlgebra:-Determinant(H71);  ## 0. - (1.89398211168871707929210978553*10^(-30))*I

Inverting the first matrix H51 results in the error. For the second we get in shortened form:

evalf[5](convert(H71^(-1),listlist));
## Result:
## [[1.1548*10^16 - (1.3826*10^16)*I, -6.7969*10^15 + (8.1376*10^15)*I], [0. - (1.8014*10^16)*I, 0.15598 + (1.0603*10^16)*I]]

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