How to test PLL for step input change in phase
I think you need to run the simulation longer, at a minimum (there could also be issues with the code).  But take for example the image I attached below.



You can see the input (top signal) and output (bottom signal).  The bottom signal does not even begin to lock on until the second period.  Even then it takes some time for it to properly lock on (10+ cycles, not shown).

I would also just implement theta as a function of time, rather than mapping it to the linear indices of time (it's simpler to me):

F=60;
step=1/(f*4)
t=0:step:1; #240 datapoints, 60 cycles of runtime
theta(t)=t*720; #linear phase shift of 2 cycles over 60 cycles
in=sin(2*pi*F*t+theta);
out=pll(in);
plot(t,out,t,theta)


I don't have access to MATLAB so I can't check if this is correct.
Simulation results would be nice to help debug as well.