As BR notes, the $q$-expansion presentation of Eisenstein series is unhelpful for determining what happens under $z\rightarrow -1/z$. (The weight-two aspect creates some additional complications, but these are not the crucial ones here. Let's ignore convergence issues, at least for a while.)

Another way to describe level-$N$ Eisenstein series is as
$$
E_{c_o,d_o,N}(z) = \sum_{c=c_o (N),\;d=d_o (N)} {1\over (cz+d)^{2k}}
$$
where $2k>2$ for convergence. It's immediate that $z\rightarrow -1/z$ interchanges $c_o$ and $d_o$, so the transformation of such Eisenstein series is clear.

For weights $>2$, the Fourier coefficients are readily computed, as in many standard sources (probably Diamond-Shurman, too). Various devices (often called "Hecke summation") extend this to the weight-two case.

Suitable linear combinations of these, summed over $c_o,d_o$ mod $N$, give the type desired. But there are complications depending on whether the characters are primitive or not mod $N$, probably a necessity of inverting matrices of Dirichlet $L$-function values, things like that. Do-able, but detail-y.

(That book-keeping is very nicely packaged by rewriting everything as modular forms on adele groups. This may be too pricey for your purposes.)

D. Bump's book on automorphic forms also does some such computations.