Eigenfilters in Matalb
It seems to me based on your question, that you are still missing some important concepts related to the algorithm, but you seem to think that the problem is not the concept. Let me see if I can clear this up.

First of all, you seem to be referencing the book Multirate Systems and Filter Banks by Vaidyanathan (pages 53-55). This would have been helpful information that you should have included in your post (or whatever other reference you obtained this from). The exposition in the book is strictly for designing eigenfilters with a low pass response (ideally with unity gain in the passband and zero gain in the stopband). The tutorial paper On The Eigenfilter Design Method and Its Applications: A Tutorial and its references extend this method to the design of filters with an arbitrary response. You should be following those references for your problem. The key point is to use reference frequencies in the frequency response to cast the cost function in the desired quadratic form so that the Rayleigh Principle can be applied to find a minimum solution (namely, the eigenvector corresponding to minimum eigenvalue).

I'd be happy to help walk you through my understanding of the filter design procedure for your particular case, but I have limited free time. For this reason, I won't do so unless I hear that you are interested.



Okay, since this is a Q&amp;A website, and I want to keep this short, I will be glossing over a lot of details that you should definitely read on your own. The tutorial paper referenced above indicates that we can obtain an arbitrary response for the eigenfilter by minimizing the cost function
$$ \xi = \mathbf{h}^\dagger \mathbf{P} \mathbf{h} $$
where I have adhered to the notation in the paper so $\dagger$ indicates the conjugate transpose (a.k.a. hermitian). Ignoring the constants and the weighting function from the paper, the matrix $\mathbf{P}$ can be computed as 
$$ \mathbf{P} = \int_{\mathcal{R}} \left[ \frac{D(\omega)}{D(\omega_0)} \mathbf{e}(e^{j\omega_0}) - \mathbf{e}(e^{j\omega}) \right]^* \left[ \frac{D(\omega)}{D(\omega_0)} \mathbf{e}(e^{j\omega_0}) - \mathbf{e}(e^{j\omega}) \right]^T d\omega $$
where $\mathcal{R}$ is the region over which $D(\omega)$ is defined, $D(\omega)$ is your desired frequency response, $\omega_0$ is a fixed reference frequency, and $\mathbf{e}(z) = \left[\begin{array}{cccc} 1 &amp; z^{-1} &amp; \cdots &amp; z^{-N} \end{array} \right]^T$. When I have more time, I will provide more details on how to go about computing this (and converting it from being complex-valued to real-valued and linear phase). Although, this is computable using trapz as you suggested in your original question, and so it may be an acceptable solution already. 



I believe that your specific case for $D(\omega)$ can be defined as follows
$$ D(\omega) = \left\{\begin{array}{cc} 1 &amp; |\omega| \leq \frac{\pi}{3}\\
0 &amp; \frac{5\pi}{12} \leq |\omega| \leq \frac{7\pi}{12} \mbox{ or } \frac{9\pi}{12} \leq |\omega| \leq \pi\\
\frac{12}{\pi}\left(|\omega| - \frac{15 \pi}{24}\right)  &amp; \frac{15\pi}{24} \leq |\omega| \leq \frac{17 \pi}{24} \\
\mbox{undefined} &amp; \mbox{otherwise} \end{array}\right.$$
where $\omega$ can vary between $-\pi$ and $\pi$. Let me know if that is correct. This gives you a total of 7 nonzero regions over which you need to integrate to obtain the elements of the matrix $\mathbf{P}$. 

A convenient reference frequency for you will be $\omega_0 = 0$ or DC. This leads to a special form where $\mathbf{e}(e^{j \omega_0}) = \mathbf{1}$ and $\mathbf{1}$ is a column vector of all ones with the same length. Also, $D(\omega_0) = 1$ so that simplifies things a bit. For now, if you choose a spectrum with conjugate symmetry, you should obtain a real response. It would be best to re-derive the formula in terms of $\mathbf{c}(\omega)$ (a vector of cosines analogous to the paper's definition for $\mathbf{e}(\omega)$ that exploits the symmetry in the coefficients). You might attempt to use this information together with the equation for the matrix $\mathbf{P}$ to find your filter coefficients. If you do and post the MATLAB that you generate, then I can do more to assist you with this.