Plot normalized phase velocity Vp/c Taflove

lunda0=1;
Dx=lunda0/5;
c=3e8;
omiga=(2*pi*c)/lunda0;
Dt=Dx/(2*c);
 alpha=1:1:90;
k=2*pi;
A=Dx*cos(1:1:90)/2;
B=Dx*sin(1:1:90)/2;
C=((Dx/(c*Dt))^2)*(sin(omiga*Dt/2))^2;
   for i=1:3
k=k-((sin(A.*k)).^2+(sin(B.*k)).^2-C)./(A.*sin(2.*A.*k)+B.*sin(2.*B.*k))

   end
kfinal(alpha)=k;
Vp=2*pi*c./kfinal(alpha);
x=Vp./c;

plot(alpha,x)