This paper presents an explicit integration scheme to compute the stiffness matrix of a nine node linear convex quadrilateral element of Lagrange family using symbolic mathematics and discretisation of polygonal domain by such finite elements using a novel auto mesh generation technique , In finite element analysis, the boundary value problems governed by second order linear partial differential equations, the element stiffness matrices are expressed as integrals of the product of global derivatives over the linear convex quadrilateral region. These matrices can be shown to depend on the material properties and the matrix of integrals with integrands as rational functions with polynomial numerator and the linear denominator (4+ ) in the bivariates over a nine node 2-square (-1 ) with nodes at the points {(-1,-1),(1,-1),(1,1),(-1,1),(0,-1),(1,0).(0,1),(- 1,0),(0,0)} in local parametric space. In this paper, we have computed these integrals in exact forms using the symbolic mathematics capabilities of MATLAB. The proposed explicit finite element integration scheme can be applied to solve boundary value problems in continuum mechanics over convex polygonal domains. We have also developed a novel auto mesh generation technique of all nine node linear convex quadrilaterals for a polygonal domain which provides the nodal coordinates and element connectivity. We have used the explicit integration scheme and this novel auto mesh generation technique to solve the Poisson equation u ,where u is unknown physical variable and in with given Dirichlet boundary conditions over the given convex polygonal domain.