Acting scalar variables on matrices while creating density plots

Hi, I am trying to make a density plot in Matlab. My problem is that I want to use some Hamiltonian,

H=a*kron(PauliX,PauliZ)+b*kron(PauliZ,PauliY)

and take the unitary evolution

U=expm(-i*H*t)

over some two qubit input state

rhoe=[0.5 0; 0 0.5]
rhon=[0.5 0; 0 0.5]
rhoin=kron(rhoe,rhon)

such that

rhoout=U*rhoin*U'

The variables in my density plot will be a, t and the value given in element 1,1 of rhoout. I want to be able to vary both a and t but I am having the problem that the matrix dimensions are wrong. Ideally, I would do a for loop for a=0:n and t=0:m but I don't know how to get this to work and also save the outputs of each rhoout. I have been trying to use the pcolor function but my problem is in getting a list of rhoout values. Can anyone tell me how I could go about coding this?