\begin{equation}\label{Parabolic} \frac{\partial \phi(\mathbf{x},t)}{\partial t} - \Delta\phi(\mathbf{x},t)+\frac{f(\phi(\mathbf{x},t))}{\epsilon^2}=0 \end{equation}
\begin{equation}\label{boundary} \frac{\partial \phi(\mathbf{x},t)}{\partial \mathbf{n}}=0, \ \text{on} \ \partial \Omega \times [0,T] \end{equation}
\begin{equation}\label{initail} \phi(\mathbf{x},0) = \phi_{0}(\mathbf{x}) \ \text{in} \ \ \Omega \end{equation} Where $f(\phi)=\phi-\phi^{3}$, $\epsilon$ is a positive constant, and $\mathbf{n}$ is the outward normal vector at the domain boundary.
I use finite element method to solve this equation:
Firstly, I discrete the equation in time space \begin{equation}\label{timediscrete} \frac{\phi^{n+1}-\phi^{n}}{\tau}-\Delta\phi^{n+1}=-\frac{f(\phi^{n})}{\epsilon^2} \end{equation} This is not a strict back Euler method, we get a elliptic equation with homogeneous neumann boundary condition： \begin{equation}\label{ellipse} \frac{\phi^{n+1}}{\tau}-\Delta\phi^{n+1}=-\frac{f(\phi^{n})}{\epsilon^2}+\frac{\phi^{n}}{\tau} \end{equation}
Suppose $V_{h}\in H^{1}(\Omega)$ is the finite element space, the problem can be stated as follows: find $\phi_{h}(\cdot,t)\in V_{h}$ s.t. \begin{equation}\label{variation} \frac{1}{\tau}(\phi_h^{n+1},v_h)-a(\phi_h^{n+1},v_h)=(g^n,v_h) \end{equation} where $v_h \in V_h, a(u,v)=\int_{\Omega}\nabla u\cdot \nabla v dxdy,\ g^n=-\frac{f(\phi^{n})}{\epsilon^2}+\frac{\phi^{n}}{\tau}$.
Assume $$\phi_{h}(\mathbf{x},t)=\sum_{i=1}^{N}\mu_{i}(t)\varphi_{i}(\mathbf{x})$$, put this in the above and set $v_h=\varphi_{i}$, we can get the following equations:
\begin{equation}\label{equationt} \sum_{i=1}^{N}(\frac{1}{\tau}(\varphi_i,\varphi_j)+a(\varphi_i,\varphi_j))\mu_i(t_{n})=(g^n,\varphi_j) \end{equation} with initial condition: \begin{equation}\label{equation0} \sum_{i=1}^{N}(\varphi_i,\varphi_j)\mu_i(0)=(\varphi_0,\varphi_j) \end{equation} $j=1,2,\cdots,N$.
As I don't know the exact solution of the problem, How can I verify my numerical method? I use $k$th degree polynomial as the basis function, is there any theoretical convergence order of my numerical method?
In order to verify your numerical method you can choose one of two things.
Instead solving $Lu=0$ solve $Lu=f$ subject to inhomogeneous boundary condition $Bu=g$. Next, choose an arbitrary function $u$ and calculate $f=Lu$ and $g=Bu$ analytically.
For example for exact solution $u=e^{at}e^{bx+dy}$ on rectangular domain you use $$f = \left(a - b^2 - d^2 +\frac{1-e^{2at}e^{2(bx+dy)}}{\epsilon^2}\right)e^{at}e^{bx+dy} $$ subject to $$ \phi_x(x_0,y,t)=be^{at}e^{bx_0+dy},\\ \phi_x(x_n,y,t)=be^{at}e^{bx_n+dy},\\ \phi_y(x,y_0,t)=de^{at}e^{bx+dy_0},\\ \phi_y(x,y_n,t)=de^{at}e^{bx+dy_n},\\ $$ and $$ \phi(x,y,0) = e^{bx+dy} \ \text{in} \ \ \Omega $$
Instead of testing against an unknown exact solution you refine the grid and measure sort of Cauchy convergence, i.e. $e=\|u_h-u_{h/2}\|$ where $u_h$ is a solution obtained on grid of size\step $h$. Measure the rate of convergence of this error which will be the same as for $e=\|u-u_h\|$ where $u$ is exact solution if known. You will need odd number of points of the grid to be able to match between subsequent grids (the finer one has twice more points, but you test only for those that a common for both grids).