Error: NaN/Inf breakpoint hit for symvar.m on line 33. 33 n = inf;

I'm trying to calculate the discrete wavenumber and plot the phase error for all the values of p from 0 to 180. I know there is a probelm with the differentiation operator but I dont know how to fix it.

My Code;

close all

syms k

c = 3*(10^8);

w = 2.0*pi*c;

ta = input('Enter the theta value ');

CFLN = input('Enter the CFLN value ');

dx = input('Enter the dx value ');

dy = input('Enter the dy value ');

dz = input('Enter the dz value ');

b = sqrt((1/dx)^2+(1/dy)^2+(1/dz)^2);

dt = CFLN /b;

for p = 0:180

sxa = (cos(p)*sin(ta)*dx)/2.0;

sya = (sin(p)*sin(ta)*dy)/2.0;

sza = (cos(ta)*dz)/2.0;

St = (1/dt)* sin((w*dt)/2.0);

sx = (1/dx)^2*sin(k*sxa)^2;

sy = (1/dy)^2*sin(k*sya)^2;

sz = (1/dz)^2*sin(k*sza)^2;

Sx = diff(sx);

Sy = diff(sy);

Sz = diff(sz);

fk = ((1/c^2)*St^2 -(sx^2+sy^2+sz^2));

fkp = -(Sx+Sy+Sz);

k = (2.0*pi) - (fk/fkp);

phase_error = -180/pi * (fk/fkp);

end

plot(phase_error,p);

My error:

Enter the theta value 90

Enter the CFLN value 0.99

Enter the dx value 0.1

Enter the dy value 0.1

Enter the dz value 0.1

NaN/Inf breakpoint hit for symvar.m on line 33.

33 n = inf;