Finding rate of convergence by curve fitting in Matlab

I have some data: number of nodes $N$ and error in energy norm corresponing to it.

I have seen in some references that the rate of convergence is reported by

$$\| u-u_h\| _E=CN^{\alpha}$$

How can I find $C$ and $\alpha$ by MATLAB, I tried by :

polyfit $(N,\| u-u_h\| _E,1)$

But the answer was not true.

How can I find them.

First of all, rates of convergence are usually given in the form $$\|u-u_h\| \leq C N^\alpha,$$ rather than equality. Furthermore, rates are asymptotic, i.e., only have to hold for $N\to \infty$. This means that you're unlikely to find a single $C$ and $\alpha$ such that your equation holds.
Another reason why your approach doesn't work is that what you're trying to fit is not a polynomial, since $\alpha$ is a) unknown and b) not an integer (for one thing, it must be negative since the error goes down as $N$ goes up).
What people usually do is look at a double-logarithmic plot: If you take the logarithm of your equation (or my inequality), you get $$log(\|u-u_h\|) \leq \log(CN^\alpha) = \log(C) + \alpha\log(N).$$ This is a linear polynomial $ax+b$ in $\log (N)$ with coefficients $a=\alpha$ and $b=\log(C)$ (from which you can find $C=\exp(b)$).
So if you apply polyfit to $\log(N)$, $\log(\|u-u_h\|)$, you should get an array $[a,b]$ with $\alpha=a$ and $C=\exp(b)$.
I think that you should be looking for a rate of convergence in the form $$\|u-u_h\|_E \leq Ch^{\alpha}$$ where $h$ is characteristic mesh size. In your case, a rough estimate would be $$h = \frac{1}{\sqrt{N}}$$ in two dimensions.
You need to evaluate the discretization error on several grids with decreasing $h$ and then use the polynomial fit as [alpha, C] = polyfit(lc, error, 1) where lc is a vector of logarithms of characteristic lengths lc = [log(h1), log(h2), log(h3), ..., log(hn)] and error is a vector of logarithms of corresponding error norms.