solving dispersion equation numerically using fsolve

I am trying to solve the dispersion equation numerically using fsolve.

M=0.4;
a=0.25;
D=2.34e6;
mu=396;
fun = @(kz) sqrt(kz.^2 - (omega / c0).^2) .* ((D * kz.^4 - mu * omega^2) .*...
            sqrt((omega / c).^2 - (1 - M^2) * kz.^2 - 2 * omega * kz / c * M) .* ...
            sin(sqrt((omega / c).^2 - (1 - M^2) * kz.^2 - 2 * omega * kz / c * M) * a) + ...
            rho * omega^2 .* (1 - kz./(omega / c) * M).^2 .* cos(sqrt((omega / c).^2 - ...
            (1 - M^2) * kz.^2 - 2 * omega * kz / c * M) * a))- rho0 * omega^2 .* ...
            sqrt((omega / c).^2 - (1 - M^2) * kz.^2 - 2 * omega * kz / c * M) .* ...
            sin(sqrt((omega / c).^2 - (1 - M^2) * kz.^2 - 2 * omega * kz / c * M) * a);

Here,M=0.4, a=0.25, D=2.34e6,mu=396. kz is the root of equation and it can be complex. I want to find the roots for general initial guesses as I dont know the asymptotic solution.for each value of omega there will be infinite values of kz as it is a transcendental equation. c-c0 and rho-rho0 are speed of sound in inner-outer medium and density of inner-outer mediums. The inner-outer mediums can be air-water, water-air,air-air,water-water. I want to plot kz vs omega for all these four combinations. The range of omega is 0 to 1e4. c_water=1500, c_air=340, rho_water=1000, rho_air=1.21.

[Torsten]

You can only find the possible values for kz if all other variables in your equation are known as numerical values.

A general symbolic solution is not possible.

[MANISH DWIVEDI]

Hey Torsten,

Thanks for responding. I need the numerical solutions only. I have mentioned all the variables in question. omega is varying from 0 to 1e4. Now for c-c0 and rho-rho0, let's say the combination of fluids (inner-outer) is air-water then c=340, c0=1500, rho=1.21 and rho0=1000. Smilarly c -c0 and rho-rho0 are defined for other cobinations of inner-outer fluids.

To be precise I want to plot Cz vs Omega for all four combinations as shown in figure 2. where Omega=omega * a/c0 and Cz =omega/kz/c0.

[Torsten]

Here is an example:

M=0.4;
a=0.25;
D=2.34e6;
mu=396;
c=340;
c0=1500;
rho=1.21;
rho0=1000;
omega =2;
fun = @(kz) sqrt(kz.^2 - (omega / c0).^2) .* ((D * kz.^4 - mu * omega^2) .*...
            sqrt((omega / c).^2 - (1 - M^2) * kz.^2 - 2 * omega * kz / c * M) .* ...
            sin(sqrt((omega / c).^2 - (1 - M^2) * kz.^2 - 2 * omega * kz / c * M) * a) + ...
            rho * omega^2 .* (1 - kz./(omega / c) * M).^2 .* cos(sqrt((omega / c).^2 - ...
            (1 - M^2) * kz.^2 - 2 * omega * kz / c * M) * a))- rho0 * omega^2 .* ...
            sqrt((omega / c).^2 - (1 - M^2) * kz.^2 - 2 * omega * kz / c * M) .* ...
            sin(sqrt((omega / c).^2 - (1 - M^2) * kz.^2 - 2 * omega * kz / c * M) * a);
kz=-1:0.001:1;
[m,idx] = min(abs(fun(kz)))

m = 0.0047

idx = 996

plot(kz,abs(fun(kz)))

I tried various values for omega, but couldn't find a parameter constellation where your function has a zero.