Clearly, for a long (L >> D) coil, this expression is valid
H = N * I
Your expression is only valid, if µ is constant over the complete region (in- and outside the coil. In practice, it only applies for µr = 1.
Because the expression doesn't contain a D parameter, it would be valid also for multiple layers of different diameter. But I doubt, that it helps with practical solenoid problems. In the general case, you can integrate Biot-Savarts rule to get the field of arbitrary rotation symmetrical coil geometries. But it still doesn't describe fields of coils with inhomogenous µr, e.g. any solenoid with a core. Except for very simple cases, numerical methods are required for an exact solution.
The right forum would be Mathematics & Physics or Electronic Elementary Questions, I think.