Solving parabolic equation using an implicit method has an advantage from a stability point of view, but on the other hand in practical engineering problems a large global system of equations must be solved in each time step. To overcome this problem the domain decomposition leads to solve instead the large original problem several smaller ones. This technique enables also to use parallel computations. This paper is a continuation of {\it C. N. Dawson} and {\it T. F. Dupont}’s explicit/implicit non-overlapping domain decomposition procedures [Math. Comput. 58, No. 197, 21--34 (1992;

Zbl 0746.65072)] and presents two kinds of new Galerkin non-overlapping domain decomposition methods which differs in a prediction of the fluxes on the boundary $\Gamma$ between subdomains. In the first method (IM-PGDD) the explicit flux calculation is just the mean integral value of the directional derivative of the solution on the boundary $\Gamma$. The second method (EIM-PGDD) is based on an extrapolation of the flux calculation in order to improve the high-order accuracy with respect to the parameter $H$ which represents the size of the boundary region $\Gamma$. Both methods are presented on a rectangular computational domain by decomposing it into two subdomains for solving a parabolic equation with a homogeneous Neumann boundary condition but the extension for other boundary conditions is included too. Then for the IM-PGDD method an error estimation of the type $$\max_{1\leq n \leq N} \| u^n -U^n \|_{L_2(\Omega)} \leq C ( \Delta t + h^{r+1} + H^{\frac52})$$ is derived, where $r$ represents polynomial degree of finite elements. For the scheme (EIM-PGDD) the final result is $$\max_{1\leq n \leq N} \|u^n -U^n \|_{L_2(\Omega)} \leq C ( \Delta t + h^{r+1} + H^{\frac92} ),$$ both under some relation on $\Delta t $ and $H$. Finally, numerical experiments that confirm the accuracy and efficiency of both presented schemes are presented.