A multidimensional generalization of Melvin’s solution for an arbitrary simple Lie algebra G is considered. The gravitational model in D dimensions, D≥ 4 , contains n 2-forms and l≥ n scalar fields, where n is the rank of G. The solution is governed by a set of n functions Hs(z) obeying n ordinary differential equations with certain boundary conditions imposed. It was conjectured earlier that these functions should be polynomials (the so-called fluxbrane polynomials). The polynomials Hs(z) , s= 1 , … , 6 , for the Lie algebra E6 are obtained and a corresponding solution for l= n= 6 is presented. The polynomials depend upon integration constants Qs, s= 1 , … , 6. They obey symmetry and duality identities. The latter ones are used in deriving asymptotic relations for solutions at large distances. The power-law asymptotic relations for E6-polynomials at large z are governed by the integer-valued matrix ν= A- 1(I+ P) , where A- 1 is the inverse Cartan matrix, I is the identity matrix and P is a permutation matrix, corresponding to a generator of the Z2-group of symmetry of the Dynkin diagram. The 2-form fluxes Φ s, s= 1 , … , 6 , are calculated. © 2017, The Author(s).