astronomy - Units after solving TOV equation - Stack Overflow

I wrote a code that solves the TOV (Tolman-Oppenheimer-Volkoff) equation. I set G=c=1, and I don't understand what units the output values have (mass and radius).

To check my code, I use a well-known solution (see image), but my problem is that I don't understand where the units of kilometers (km) come from, especially in relation to my code.

The graph of the known solution looks exactly the same as mine. The differences are that r is given in km and M in units of M_star (solar mass). However, I don't understand where in the code I specified that r should be output in km.

#include <stdio.h>
#include <math.h>
#include <gsl/gsl_odeiv2.h>
#include <gsl/gsl_errno.h>

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

double eos(double p)
{
    if (p < 1e-20) return 0.0; 
    //if (p > 1e20) return 1e20; // limit pressure
      return pow(p / 20.0, 3.0 / 5);
}

int func(double r, const double y[], double f[], void *params)
{

    double e = eos(y[0]);



    // Normale TOV-Gleichungen
        if (fabs(1 - 2 * y[1] / r) < 1e-13)
        {
        return GSL_FAILURE; // Avoid Singularity
        };
        

    f[0] = -(e + y[0]) * (y[1] + 4 * M_PI * pow(r, 3) * y[0]) / (r * (r - 2 * y[1])); // dp/dr = ...
    f[1] = 4 * M_PI * r * r * e;  // dm/dr=...

    return GSL_SUCCESS;

}

int main()
{




    FILE *file = fopen("output.txt", "w");
    if (file == NULL)
    {
        printf("Fehler beim Öffnen der Datei!\n");
        return 1;
    }

    gsl_odeiv2_system sys = {func, NULL, 2, NULL};
    gsl_odeiv2_driver *driver = gsl_odeiv2_driver_alloc_y_new(
        &sys, gsl_odeiv2_step_rk4, 1e-15, 1e-15, 1e-15);

     double r = 1e-6;         // Starting radius
    double y[2] = {1e-13, 0}; // Initial conditions: {pressure, mass}

    // fprintf(file, " r       P       m\n");

    for (int i = 0; i <= 130; i++)
    {

        r = 1e-6;
        y[0] = 1e-13 * pow(1.8, i); 
        y[1] = 0;                   

        fprintf(file, " P_0 = %e ", y[0]);


        while (y[0] > 1e-17)
        {
            // fprintf(file, "%lf %e %lf\n", r, y[0], y[1]);
            int status = gsl_odeiv2_driver_apply(driver, &r, r + 0.001, y);
            if (status != GSL_SUCCESS)
            {
                fprintf(stderr, "Error: %s\n", gsl_strerror(status));
                break;
            }
        }

        fprintf(file, " r = %e M= %e\n", r, y[1]); 
        gsl_odeiv2_driver_reset(driver); // Reset driver for each iteration

    }

    gsl_odeiv2_driver_free(driver);
    fclose(file);

    return 0;
}

I've looked at the code many times, but I can't figure it out. I'm using geometric units (with G=c=1), but when I try to apply them to my radius (convert them to kilometers), I suddenly get different results compared to the known solution.