A finite element formulation for bending analysis of isotropic and orthotropic plates based on two-variable refined plate theory is developed in this paper. The two-variable refined plate theory which can be used for both thin and thick plates predicts parabolic variation of transverse shear stresses across the plate thickness and therefore, it does not need shear correction factor in the formulation and the zero stress conditions are satisfied on free surfaces. The von-Karman nonlinear terms are considered in strain-displacement equations and governing equations are derived using the Hamilton's principle. After constructing weak form equations, a new 4-node rectangular plate element with six degrees of freedom at each node is used for discretization of the domain. The non-linear coupled governing equations are solved by Newton–Raphson method. The finite element code is written in MATLAB which can be used for analysis of thin and thick, isotropic and orthotropic plates with various boundary conditions. Some benchmark problems are solved by the developed code and the obtained displacements and stresses are compared with the existing results in the literature which show the accuracy and efficiency of presented finite element formulation.