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.

Make sense?