I have written a physics simulation about the classical mechanics famous problem, the three body problem. I have written some code in c++ and have constructed a class named ‘System_variables’.

#include <math.h>
#include <iostream>
using namespace std;

class System_variables
{
    private:

        // These are the system variables:- position, velocity and accleration
        long double body_1_pos[3];
        long double body_2_pos[3];
        long double body_3_pos[3];

        long double body_1_v[3];
        long double body_2_v[3];
        long double body_3_v[3];

        long double acc_body_1[3];
        long double acc_body_2[3];
        long double acc_body_3[3];

        // Needed constants
        const double body_mass = 6 * pow(10, 24);
        double time_step;
        const float G = 6.67 * pow(10, -11);

    public:


        System_variables(long double body_1_pos_i[3], long double body_2_pos_i[3], long double body_3_pos_i[3], long double body_1_v_i[3], long double body_2_v_i[3], long double body_3_v_i[3])
        {
            // Initiates the system variables, position and velocity.
            for (int m = 0; m>=2; m++)
            {
                body_1_pos[m] = body_1_pos_i[m];
                body_2_pos[m] = body_2_pos_i[m];
                body_3_pos[m] = body_3_pos_i[m];

                body_1_v[m] = body_1_v_i[m];
                body_2_v[m] = body_2_v_i[m];
                body_3_v[m] = body_3_v_i[m];
            }
        }
        void init_time_step(double step)
        {
            // Initializes time step.
            time_step = step;
        }
        void calculate_acc()
        {   
            // Calculates the accleration on the bodies using newtons equations.
            for (int m=0; m<=2; m++)
            {
                acc_body_1[m] = G * body_mass/pow((body_2_pos[m] - body_1_pos[m]), 2) + G * body_mass/pow((body_3_pos[m] - body_1_pos[m]), 2);
                acc_body_2[m] = G * body_mass/pow((body_1_pos[m] - body_2_pos[m]), 2) + G * body_mass/pow((body_3_pos[m] - body_2_pos[m]), 2);
                acc_body_3[m] = G * body_mass/pow((body_1_pos[m] - body_3_pos[m]), 2) + G * body_mass/pow((body_2_pos[m] - body_3_pos[m]), 2);
            };
            return;
        }

        void apply_acc_to_v()
        {   
            // Applys the accleration to the velocities of the bodies.
            for (int m=0; m<=2; m++)
            {
                body_1_v[m] += acc_body_1[m] * time_step;
                body_2_v[m] += acc_body_2[m] * time_step;
                body_3_v[m] += acc_body_3[m] * time_step;
            }
            return;
        }

        void apply_v_to_pos()
        {   
            // Applys the velocities to the positions of the bodies.
            for (int m=0; m<=2; m++)
            {
                body_1_pos[m] += body_1_v[m] * time_step;
                body_2_pos[m] += body_2_v[m] * time_step;
                body_3_pos[m] += body_3_v[m] * time_step;
            }
            return;
        }

        long double get_body_1_pos(int index)
        {
            // Returns the position of body 1 based on index.
            return body_1_pos[index];
        }

        long double get_body_2_pos(int index)
        {   
            // Returns the position of body 2 based on index.
            return body_2_pos[index];
        }

        long double get_body_3_pos(int index)
        {
            // Returns the position of body 3 based on index.
            return body_3_pos[index];
        }

        long double get_body_1_v(int index)
        {
            // Returns the velocity of body 1 based on index.
            return body_1_v[index];
        }

        long double get_body_2_v(int index)
        {
            // Returns the velocity of body 2 based on index.
            return body_2_v[index];
        }

        long double get_body_3_v(int index)
        {
            // Returns the velocity of body 3 based on index.
            return body_3_v[index];
        }


};

int main()
{
    long double array_1[3] = {100,100,100}; 
    long double array_2[3] ={10000000000, 10000000000, 10000000000};
    long double array_3[3] = {10000000000, 10000000000, -10000000000};
    long double array_0[3] =  {0,0,0};

    System_variables system(array_1, array_2, array_3, array_0, array_0, array_0);
    system.init_time_step(1);

    for (int time = 8; time <= 10; time++)
    {
        cout << "x position of body 1: " << system.get_body_1_pos(0) << endl 
        << "y position of body 1: " << system.get_body_1_pos(1) << endl 
        << "z position of body 1: " << system.get_body_1_pos(2) << endl << endl;

        cout << "x position of body 2: " << system.get_body_2_pos(0) << endl 
        << "y position of body 2: " << system.get_body_2_pos(1) << endl 
        << "z position of body 2: " << system.get_body_2_pos(2) << endl << endl;

        cout << "x position of body 3: " << system.get_body_3_pos(0) << endl 
        << "y position of body 3: " << system.get_body_3_pos(1) << endl 
        << "z position of body 3: " << system.get_body_3_pos(2) << endl << endl;

        system.calculate_acc();
        system.apply_acc_to_v();
        system.apply_v_to_pos();
    }
}

I thought I did everything fine but when I ran the program…

x position of body 1: 8.60984e-305
y position of body 1: 3.39519e-313
z position of body 1: 9.34072e-305

x position of body 2: 1.45708e-303
y position of body 2: 1.69826e-313
z position of body 2: 6.92657e-303

x position of body 3: 9.96985e-317
y position of body 3: 2.09368e+268
z position of body 3: 0

x position of body 1: inf
y position of body 1: inf
z position of body 1: inf

x position of body 2: inf
y position of body 2: inf
z position of body 2: inf

x position of body 3: inf
y position of body 3: 2.09368e+268
z position of body 3: inf

x position of body 1: -nan(ind)
y position of body 1: -nan(ind)
z position of body 1: -nan(ind)

x position of body 2: -nan(ind)
y position of body 2: -nan(ind)
z position of body 2: -nan(ind)

x position of body 3: -nan(ind)
y position of body 3: 2.09368e+268
z position of body 3: -nan(ind)

It returned the following. Why is it outputting so much jibbrish even in the first time step? I have used the long double container, even though the numbers break once they get smaller than 3.4 * 10^(-34). How can I solve this?