The formulas you're seeing are for the infinite solenoid, which doesn't help much.
Try:
https://paginas.fe.up.pt/~ee08173/wp-content/uploads/2014/03/finite-solenoid.pdfLooks like the axial solution avoids scary integrals and simplifies to an Excel-able formula. Note, the brackets-sub/superscript notation I believe is saying you need to take the difference of the expression evaluated at the two values of ξ (and ξ is given elsewhere).
For a flat field over a region of space (i.e., dB/dx --> 0, and higher derivatives and other axes if possible), consider Helmholtz coils. There are at least two well-known cases (2 and 4 loops, giving 1st and 2nd derivatives, respectively). This will be much flatter than a solenoid.
You could also -- since this is a symbolic solution, take the derivative, and put two solenoids together such that their derivatives exactly cancel at the point of interest. Making a Helmholtz coil with finite length loops. May allow higher flux density for a given amount of wire and power dissipation?
I suppose ideally you'd want to integrate over a range of
r as well, to get a range of radii, modeling a finite thickness solenoid -- a multilayer one, perhaps.
But yeah, things quickly get complicated in that direction, and you'll find it very attractive to get more tools, or simply model everything in the first place. You might pick up Maxima to automate symbolic math. You can perform a numerical integration in Excel by entering the formula across a large enough array of cells, and summing their results -- Riemann summation. Or, learn FEMM, as mentioned.

Tim