Based on the equation you posted, the variable A=2*pi needs to be squared when it's inside the square root. But your code doesn't square A just f.
It's a messy formula so you should double check a couple of points by hand and then check for example the values
symmetric([0 3.5],[0 22])
Also you can put it in a separate function to avoid using 2 functions
fimplicit(@symmetric,[0 3.5 0 22])
function out = symmetric(cp,f)
CL = 5.850;% Longitudinal velocity in [mm/µseg]
CT = 3.184;% Shear velocity in [mm/µseg]
h=1;
om = 2*pi.*f;
k = om./cp;
p = sqrt(om.^2/CL^2-k.^2);
q = sqrt(om.^2/CT^2-k.^2);
out = tan(q*h)./q + 4*k.^2.*p.*tan(p*h)./(q.^2-k.^2).^2;
end
I noticed the plot changed a little when I use this method however.
