I have tested this model with RK method of order 8 to understand how it depends on input error and step size as well. First we should note, that relativistic model has hidden suggestion in the form $c=1$ where c is the speed of light. In this case the unit of time and coordinate is the same, for instance, meter. Second, we take Mercury semimajor axis, perihelion, eccentricity, period in days and maximal orbital speed as
Mmercury = 3.3010 10^23; RAmercury = 5.7909227 10^10; RPmercury = 4.6002 10^10; Emercury = 0.20563593; Pmercury = 87.97; vp=58983;
We use Mercury perihelion as a scale of length, therefore in this case we have parameter M of the model and initial condition as follows
c0 = 299792458.; GM = 1.32712 10^(20); mS = GM/c0^2; M=mS/RPmercury ;
$$t(0)= 0, r(0) = 1, \phi(0) = 0, t'(0) = u_1,
r'(0) = 0, \phi'(0)= vp/c0$$
Here u1 is computed as a root of equation
$$u.g.u=0, u=(u1,0,0,vp/c0)$$
g is the metric tensor. We have calculated, that $u1=1.0000000514518859$. Also we define time scale $t_0=RPmercury/c0=153.446 s$, hence in equations we should normalize $t$ on $t_0$. Finally we compute orbit as shown in Figure 1
According this calculation the period is about $49569*t0=7.60617*10^6$, while astronomical data says that period is about $Pmercury *24*3600=7.60061*10^6$. Therefore we have discrepancies in time period, but it is not due to relativistic dilation , but due to error in definition of input parameter and computation method.