The Symbolic Math Toolbox may not be appropriate for this problem.
Try this:
freq = linspace(10*1e3, 31*1e6,100);
w = 2*pi.*freq; % angular frequency
p= 6.5 *1e-3;
d= 0.6*1e-3;
r = d/2;
mu0 = 4*pi*1e-7;
mur= 200;
sigw= 8.2183e+05;
Rw = (pi*sigw*(r^2))^(-1);
tau = mu0*mur*sigw*(r^2);
num = sqrt(1i.*w.*tau).*besseli(0, sqrt(1i.*w.*tau));
denum = 2.*besseli(1, sqrt(1i.*w.*tau));
Zw = Rw.*num./denum;
Ls = -((mu0*(p+d))/(2*pi)).*(log(1-exp(-(pi*d)./(d+p)))); %^(-1);
Zs = Zw*(d+p)+1i.*w.*Ls;
% syms sigma
for k = 1:numel(w)
Zs2 = @(sigma) (1+1j) .* ((w(k).*mur.*mu0) ./ (2.*sigma)).^(1/2) - Zs(k);
sigmav(k) = fsolve(Zs2, rand*1+1i);
end
figure
plot(w, real(sigmav))
hold on
plot(w, imag(sigmav))
hold off
grid
xlabel('\omega')
legend('\Re \sigma', '\Im \sigma', 'Location','best')
Experiment to get different results.
