After a long searching for the analytic form of finite solenoid field, with the help of  from :

  1. Edmund E. Callaghan, S. H. (1960). The Magnetic Field of a Finite Solenoid (Techical note D-465). Washington, USA: Nation Aeronautics and Space Administration.
  2. Jackson, J. D. (1998). Classical Electrodynamics. John Wiley & Sons, Inc.
  3. Milton Abramowitz, I. A. (1965). Handbook of mathematical functions : with formulas, graphs, and mathematical tables. Dover.
  4. NIST Digital Library of Mathematical Functions. (n.d.). Retrieved from http://dlmf.nist.gov/
we have the analytic form.
the potential is:
A_\phi = \frac{\mu_0 I}{2\pi } \frac{1}{L} \sqrt{\frac{a}{\rho}} \left[ \zeta k \left( \frac{k^2+h^2-h^2k^2}{h^2k^2}K(k^2)-\frac{1}{k^2}E(k^2) +\frac{h^2-1}{h^2} \Pi(h^2,k^2) \right) \right]_{\zeta_-}^{\zeta_+}
the field vector is:
\displaystyle B_\rho = \frac{\mu_0 I}{2\pi} \frac{1}{L} \sqrt{\frac{a}{\rho}} \left[ \frac{k^2-2}{k}K(k^2) + \frac{2}{k} E(k^2)\right]_{\zeta_-}^{\zeta_+}
\displaystyle B_z =\frac{\mu_0 I}{2\pi} \frac{1}{L} \frac{1}{2 \sqrt{a \rho}} \left[ \zeta k \left(K(k^2) + \frac{a-\rho}{a+\rho} \Pi(h^2,k^2)\right)\right]_{\zeta_-}^{\zeta_+}
where
\displaystyle \zeta_{\pm}=z\pm \frac{L}{2}
\displaystyle h^2=\frac{4a\rho}{(a+\rho)^2}
\displaystyle k^2=\frac{4a\rho}{(a+\rho)^2+\zeta^2}
and
\displaystyle K(m)=\int_0^{\pi/2}{\frac{1}{\sqrt{1-m \sin^2 \theta }}} d\theta
\displaystyle E(m)=\int_0^{\pi/2}{\sqrt{1-m \sin^2 \theta} } d\theta
\displaystyle \Pi(n,m)=\int_0^{\pi/2}{\frac{1}{(1-n \sin^2 \theta)\sqrt{1-m \sin^2 \theta }}} d\theta
those are the complete elliptic integral of 1st , 2nd and third kind. I quoted them in here because the argument is a bit confusing across different references and I settle it down on the context of this page. The below is the plotted field line and intensity for the coil dimension is 1 radius and 1 length.