Computing the RPY tensor with C++/Eigen (linear algebra)
Using p + 1 instead of just p

You declare

const int p = 4000;


but use p+1 at all the places that matter. I think you should just use

const int p = 4001;


and replace all usages of p + 1 by just p.

Simplify a few expressions

You have repeated usage of 3 * ( i + 1) - 3 and 3 * ( j + 1) - 3. They can simply be 3*i and 3*j.

Give some computations a name

You have the following expressions:

0.75 * distance_ratio * (1.0 - 2.0*distance_ratio * distance_ratio) * rvec * rvec.transpose()

0.75 * distance_ratio * (1.0 + 2.0 / 3.0 * distance_ratio * distance_ratio)

3.0 / 32.0 / distance_ratio * rvec*rvec.transpose()

1.0 - 9.0 / 32.0 / distance_ratio


I don't know that they compute but they compute four different things. It will be easier to understand your code if you can give those computations a name and put them in a function.

double firstTerm(double distance_ratio,
                 Vector3d const&amp; rvec)
{
   return (0.75 * distance_ratio * (1.0 - 2.0*distance_ratio * distance_ratio) * rvec * rvec.transpose());
}


and use it

   rpy.block&lt;3, 3&gt;(3*i, 3*j) = firstTerm(distanc_ratio, recv);


PS I am not suggesting firstTerm as a function name. It just illustrates the concept.

Avoid multiple statements in one line

The line

    coords[i](0) = dis(gen); coords[i](1) = dis(gen); coords[i](2) = dis(gen);


is less readable than:

    coords[i](0) = dis(gen);
    coords[i](1) = dis(gen);
    coords[i](2) = dis(gen);


You could also put them in a for loop to avoid code repetition.

 for ( int j = 0; i &lt; 3; ++j )
    coords[i](j) = dis(gen);


Using const int as an argument type`

I am not a fan of using const int as an argument type. Just int works without any problem. YMMV.

The sqrt line

        Rij = sqrt(rvec(0)*rvec(0) + rvec(1)*rvec(1) + rvec(2)*rvec(2));


can be modified to

double square(double in)
{
    return in*in;
}

Rij = sqrt(square(rvec(0)) + square(rvec(1)) + square(rvec(2));


I think this not only avoids code repetition but also uses the right abstraction in the computation.