Do you have a source for the equations you're trying to implement to solve for u?
I thought that the controllability Grammian as the solution of the discrecte Lyapunov equation is the sum over n = 0 to inf, as discussed in
But in this problem we only need the "portion" of the Grammian from n = 0 to 50. So the first thing is to find the corrrect equation that solves for u. Then we can figure out the Matlab code. Is this link helpful: General State Transfer and Reachabilty ?
Also, the lyap() function is for continuous systems. The dlyap() function is used for discrete systems. However, neither of these functions need be used for this problem, I don't believe.