I'm puzzled by your "The initial condition is u(0,t) = cos^2(pi x)". You've plugged in zero for z on the LHS, but not on the RHS. Did you mean "The initial condition is u(x,0) = cos^2(pi x)."? If so, the usual procedure for applying the orthogonality condition is to set your initial condition equal to the solution sum, with t=0 (assuming your work up to now is correct):
(note the required dependence on of the coefficients ), and then multiply through by an orthogonal term:
You then integrate both sides w.r.t. from to :
You exploit the orthogonality of the cosine functions to simplify the RHS. See what you get when you do all that.
This is standard Fourier analysis, by the way.