Не могу найти ошибку

#include <iostream>
#include <cmath>
#include <iomanip>
using namespace std;
 
int main(){
    setlocale(LC_CTYPE,"rus");
    double A = 5 * 1e9, n = 4, gamma = 3;
    double p_no = 7500;
    double L_no = 0.01;
    double D = 8000;
    double p_vv = 1600;
    long double R_0;
    double M_dv = 40;
    int vid_dv = 235;
    double R_no = 0.19;
    double Pi = 3.14 ;
    double r_vo;
    double r_no;
    double R_vo = R_no - L_no;
    long double U_n;
    double p_dv = 18700;
    double x = 0.05, y = 0.01, r_sh = pow(3 * M_dv / (4 * Pi * p_dv) , 1./3.);
    double r_n = r_sh + 1e-6,Ot = R_vo - x,Do = r_sh + y;
    double shag = (Do - Ot) / 5;
    double mass[6];
    int j=0; 
    cout << " ðàäèóñ øàðà =" << r_sh <<endl;
    cout << " Îò =" << Ot <<endl;
    cout << " Äî =" << Do <<endl;
    cout << " øàã = " <<shag<<endl;
    //Òàáëèöà çíà÷åíèé r_no.
    cout << "   r_no \n" << endl;
    r_no = Ot + shag;
    mass[0] = r_no;
    cout << "r_no (1) = "<< mass[0] << endl;    
    for ( int i = 1; i <= 5 ; i++){
        r_no += shag; 
        mass[i] = r_no;
        cout << "r_no ("<< i << ") = "<< mass[i] << endl;   
    }
    cout << " \n" << endl;
    for( int i = 0; i<= 5; i++){
    cout << "òî÷êà ¹" << i <<"\n"<<endl;
    // 1) Ðàññ÷èòûâàåì âåëè÷èíó R_0.
    double M_no = (4 / 3) * Pi * p_no * ( pow(R_no + L_no, 3) - pow(R_no,3));
    double M_vv = (4 / 3) * Pi * p_vv *( pow(R_no,3) - pow(mass[i],3));// Ìàññà ÂÂ.
    R_0 = (mass[i] * M_dv + Pi * p_vv * ( pow( R_no,4) - pow( mass[i],4)) + R_no * M_no) / (M_dv + M_no + M_vv);
    
    cout << " Ro = " << R_0 <<endl;
    // 2) Ðàà÷èòàåì U_n.
    double P_go = p_vv * pow(D,2) / (2 * (gamma + 1));
    double V_0 = (4 / 3) * Pi * (pow(R_0,3) - pow(mass[i],3));
    double V_ot_r_n = (4 / 3) * Pi *(pow(R_0,3) - pow(r_n,3));
    double r_v = pow( pow(r_n,3) - 3 * M_dv / (4 * Pi * p_dv) , 1./3.);
    double y_1 = y - 1;
    U_n = pow( (P_go * V_0 / (gamma - 1) * ( 1 - pow(V_0/V_ot_r_n, y_1)) / ((2 * Pi * p_dv * pow(r_n,3) * (r_n - r_v) / r_v) + ((2/5) * Pi * p_vv * V_0 / V_ot_r_n) * (pow (R_0,3)/6 + pow(R_0,2) * r_n/3 + R_0*pow(r_n,2)/2 -pow(r_n,3)))) , 1./2.);
    cout << " Uí = " << U_n <<endl;
    // 3) Ðàññ÷èòàåì Åêâä(r_í):
    double E_kvd_ot_r_n = 2* Pi * p_dv * pow(U_n,2) * pow(r_n,3) * (r_n - r_v) / r_v;
    cout << " Eêâä( r_í ) = " << E_kvd_ot_r_n <<endl;
    cout << " \n" << endl;
    }
    
    return 0;
}

#include <iostream>
#include <cmath>
#include <iomanip>
using namespace std;
 
int main(){
    setlocale(LC_CTYPE,"rus");
    double A = 5 * 1e9, n = 4, gamma = 3;
    double p_no = 7500;
    double L_no = 0.01;
    double D = 8000;
    double p_vv = 1600;
    long double R_0;
    double M_dv = 40;
    int vid_dv = 235;
    double R_no = 0.19;
    double Pi = 3.14 ;
    double r_vo;
    double r_no;
    double R_vo = R_no - L_no;
    long double U_n;
    double p_dv = 18700;
    double x = 0.05, y = 0.01, r_sh = pow(3 * M_dv / (4 * Pi * p_dv) , 1./3.);
    double r_n = r_sh + 1e-6,Ot = R_vo - x,Do = r_sh + y;
    double shag = (Do - Ot) / 5;
    double mass[6];
    int j=0; 
    cout << " радиус шара =" << r_sh <<endl;
    cout << " От =" << Ot <<endl;
    cout << " До =" << Do <<endl;
    cout << " шаг = " <<shag<<endl;
   //Таблица значений r_no.
    cout << "   r_no \n" << endl;
    r_no = Ot + shag;
    mass[0] = r_no;
    cout << "r_no (1) = "<< mass[0] << endl;    
    for ( int i = 1; i <= 5 ; i++){
        r_no += shag; 
        mass[i] = r_no;
        cout << "r_no ("<< i << ") = "<< mass[i] << endl;   
    }
    cout << " \n" << endl;
    for( int i = 0; i<= 5; i++){
    cout << "точка №" << i <<"\n"<<endl;
    // 1) Рассчитываем величину R_0.
    double M_no = (4 / 3) * Pi * p_no * ( pow(R_no + L_no, 3) - pow(R_no,3));
    double M_vv = (4 / 3) * Pi * p_vv *( pow(R_no,3) - pow(mass[i],3));
    R_0 = (mass[i] * M_dv + Pi * p_vv * ( pow( R_no,4) - pow( mass[i],4)) + R_no * M_no) / (M_dv + M_no + M_vv);
    
    cout << " Ro = " << R_0 <<endl;
    // 2) Расчитаем U_n.
    double P_go = p_vv * pow(D,2) / (2 * (gamma + 1));
    double V_0 = (4 / 3) * Pi * (pow(R_0,3) - pow(mass[i],3));
    double V_ot_r_n = (4 / 3) * Pi *(pow(R_0,3) - pow(r_n,3));
    double r_v = pow( pow(r_n,3) - 3 * M_dv / (4 * Pi * p_dv) , 1./3.);
    double y_1 = y - 1;
    U_n = pow( (P_go * V_0 / (gamma - 1) * ( 1 - pow(V_0/V_ot_r_n, y_1)) / ((2 * Pi * p_dv * pow(r_n,3) * (r_n - r_v) / r_v) + ((2/5) * Pi * p_vv * V_0 / V_ot_r_n) * (pow (R_0,3)/6 + pow(R_0,2) * r_n/3 + R_0*pow(r_n,2)/2 -pow(r_n,3)))) , 1./2.);
    cout << "Uн = " << U_n <<endl;
    // 3) Рассчитаем Еквд(r_н):
    double E_kvd_ot_r_n = 2* Pi * p_dv * pow(U_n,2) * pow(r_n,3) * (r_n - r_v) / r_v;
    cout << " Eквд( r_н ) = " << E_kvd_ot_r_n <<endl;
    cout << " \n" << endl;
    }
    
    return 0;
}