Also, instead of using Newtonian acceleration, I tried out simply using the derivative of the time dilation in the Schwarzschild solution. What I mean is that it's more natural, and it works quite fine:
void proceed_Euler(custom_math::vector_3& pos, custom_math::vector_3& vel, const long double G, const long double dt)
{
custom_math::vector_3 grav_dir = sun_pos - pos;
const long double distance = grav_dir.length();
const long double Rs = 2 * grav_constant * sun_mass / (speed_of_light * speed_of_light);
const float beta = sqrt(1.0 - Rs / distance);
const long double gradient = (speed_of_light * speed_of_light * Rs) / (2 * distance * distance * sqrt(1 - Rs / distance));
const float alpha = 2.0 - sqrt(1.0 - (vel.length() * vel.length()) / (speed_of_light * speed_of_light));
grav_dir.normalize();
custom_math::vector_3 accel = grav_dir * gradient;
vel += accel * alpha * dt;
if (vel.length() > speed_of_light)
{
vel.normalize();
vel *= speed_of_light;
}
pos += vel * beta * dt;
}