FIR filter design: evaluating the error function
There are many problems in your question and in your code. First, Eqs $(2)$ and $(3)$ are wrong. Eq. $(2)$ cannot equal Eq. $(1)$ because it is its complex conjugate (and $\mathbf{c}$ is a complex vector). Eq. $(3)$ is meaningless.

I will not rewrite your Matlab code, but you should have added that it doesn't run.
There are many errors, among which several incompatible matrix dimensions. You should use the function cfirls instead of just copying some of its code. Then vector dimensions would be checked and care would be taken that all vectors have the same format (I chose columns).

Assuming that all these problems are resolved (simply by calling cfirls with the vectors f, d, and w), you can compute the error, but this is where the next problem occurs. Why don't you just use Eq. $(1)$ in your question to compute the frequency response? (You might as well use the Matlab command freqz). The way you define the matrix c in order to compute the frequency response H doesn't make sense.

Assuming that you manage to compute the designed frequency response on a grid, you could choose to compute the unweighted error, as suggested by your code (even though you minimized a weighted error). However, your variable S will be a scalar (as soon as the code to compute it is fixed). So how do you want to plot it as a function of frequency (as suggested by the last line of your code)?

In sum, apart from using cfirls, which will give you the correct filter coefficients, think about how you should compute the frequency response, and why you defined the matrix c the way you did. Then think about which quantity you would like to plot. As mentioned before, the squared error is just a number, not a function of frequency. It could be useful to plot the magnitude of the difference between desired and designed frequency response:

$$\big|E(\omega)\big|=\big|H(\omega)-D(\omega)\big|\tag{1}$$

With N, f, d, and w defined as in your code, the most straightforward way to do what you're trying to do would be:


h = cfirls(N,f,d,w);
H = freqz(h,1,pi*f);
E = H - d;
plot(f,20*log10(abs(E)))