\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?