Abstract: Starting from the analytical solution for the non-axisymmetric elastic problem of multi-layered soils, the analytical layer element of a single soil layer, namely a symmetric exact stiffness matrix, is acquired in the Hankel transform domain. The global stiffness matrix is obtained by assembling the interrelated layer-elements based on the principles of finite layer method. By solving algebraic equations of the global stiffness matrix which satisfy the boundary conditions, we can get the solution to the non-axisymmetric loading problem for multi-layered soils in the Hankel transform domain. Eventually the actual solutions in the physical domain can be acquired by inverting the Hankel transform. Numerical calculation and analysis are carried out by means of the corresponding computer program, and the surface horizontal displacement of multi-layered soils along the radial direction under non-axisymmetric loading is discussed.