plot issue 'input function must return single or double value. sym found.'

I want to plot NLI, but error message 'input function must return single or double value. sym found.' displayed. What should I do to plot NLI?

I tried use fplot instead of plot, but same message displayed.

clear all;
a = 0.2;      %field loss coefficient α[db/km]
alfa = a*log(20)/20;
gamma = 1.3;       % fiber non-linearity coefficient　γ[1/w/km]
Ns = 20;        %number of span
Ls = 100;       %span length[km]         
beta2 = 20.7;     %dispersion coefficient[ps^2/km]
roll = 0.3; %roll-off of gwdm
%% formula of　ρ
syms f f1 f2 real %f,f1,f2 are THz
Le = (1-exp(-2*alfa*Ls))/(2*alfa);
x1(f1,f2,f) = 1-exp(-2*alfa*Ls)*exp(4j*(pi^2)*(f1-f)*(f2-f)*beta2*Ls);%symfun
x2(f1,f2,f) = 2*alfa-(4j*(pi^2)*(f1-f)*(f2-f)*beta2);
p(f1,f2,f) = Le^(-2)*(abs(x1(f1,f2,f)/x2(f1,f2,f)))^2;
%% formula of Gwdm
syms t w k 
T=1/(32e-3);
A = pi*t/T;
x(t)= (sin(A)/A)*(cos(roll*A)/(1-(2*roll*t/T)^2));
X(w) = simplify(fourier(rewrite(x(t),'exp'),t,w));
X(f) =X(2*sym(pi)*f);
X(f)= rewrite(abs(X(f)),'sqrt');
xfunc = matlabFunction(X(f));%@(f)
GWDM(f)= (symsum(xfunc(f+(50e-3)*k),k,-5,5)/((32/8.8)*T));
func = p*GWDM(f1)*GWDM(f2)*GWDM(f1+f2-f);
ft=matlabFunction(func);
%% formula of GNLI
pint = @(f)integral2(@(f1,f2)ft(f1,f2,f),-271e-3,271e-3,-271e-3,271e-3,'AbsTol',1e-6);
NLI = @(f)(16/27).*gamma.^2.*Le.^2.*pint(f); % NLI of 1span
%fplot(NLI(f), [-0.4 0.4])
plot(f,NLI(f));

if the plot success correctly, the shape of graph is near to blue line(vertical line indicate how much larger is GNLI than GWDM.)

Answer 1

[Torsten]

a = 0.2;      %field loss coefficient α[db/km]
alfa = a*log(20)/20;
gamma = 1.3;       % fiber non-linearity coefficient　γ[1/w/km]
Ns = 20;        %number of span
Ls = 100;       %span length[km]         
beta2 = 20.7;     %dispersion coefficient[ps^2/km]
roll = 0.3; %roll-off of gwdm
%% formula of　ρ
syms f f1 f2 real %f,f1,f2 are THz
Le = (1-exp(-2*alfa*Ls))/(2*alfa);
x1(f1,f2,f) = 1-exp(-2*alfa*Ls)*exp(4j*(pi^2)*(f1-f)*(f2-f)*beta2*Ls);%symfun
x2(f1,f2,f) = 2*alfa-(4j*(pi^2)*(f1-f)*(f2-f)*beta2);
p(f1,f2,f) = Le^(-2)*(abs(x1(f1,f2,f)/x2(f1,f2,f)))^2;
%% formula of Gwdm
syms t w k 
T=1/(32e-3);
A = pi*t/T;
x(t)= (sin(A)/A)*(cos(roll*A)/(1-(2*roll*t/T)^2));
X(w) = simplify(fourier(rewrite(x(t),'exp'),t,w));
X(f) =X(2*sym(pi)*f);
X(f)= rewrite(abs(X(f)),'sqrt');
xfunc = matlabFunction(X(f));%@(f)
GWDM(f)= (symsum(xfunc(f+(50e-3)*k),k,-5,5)/((32/8.8)*T));
func = p*GWDM(f1)*GWDM(f2)*GWDM(f1+f2-f);
ft=matlabFunction(func);
%% formula of GNLI
pint = @(f)integral2(@(f1,f2)ft(f1,f2,f),-271e-3,271e-3,-271e-3,271e-3,'AbsTol',1e-6);
NLI = @(f)(16/27).*gamma.^2.*Le.^2.*pint(f); % NLI of 1span
f = -0.4:0.01:0.4;
plot(f,arrayfun(@(f)NLI(f),f));

[柊介 小山内]

Thanks! I can plot graph.