Help me with plotting, please.

Here is the code:

h_s=0.1; % step for integrals and graphs
n_0=1*10^19; %quantity of atoms
W=0.7; %wide (cm)
H=0.3; %height (cm)
L=0.78; %length (cm)
alpha=5; %absorption coefficient (1/cm)
I_s=3.5; %signal intensity (W)
P_in=5; %input power (W)
sigma_a=5*10^(-19); %absorption cross-section (cm^2)
sigma_e=1*10^(-18); %emission cross-section (cm^2)
h=6.62*10^(-27); %Plank's constant (cm^2*g/s)
tau=4*10^(-6); %lifetime of the exited state (s)
lambda_0=1.9*10^-4;%wavelength of the pump (cm)
c=3*10^10; %light velocity (cm/s)
w_s=80*1^-4; %spot radius (cm)
I_sat=h*c*tau/(lambda_0*(sigma_e+sigma_a));
syms x y
P1=log(y/P_in);
P2=(sigma_e/sigma_a)*(log(1+x*sigma_a*tau*lambda_0/(h.*c))-log(1+(x*sigma_a*tau*lambda_0/(h.*c))*exp(-alpha*L)));
P3=(P_in*(pi/2)*w_s/I_sat)*(y/P_in-1);
f(x,y)=P1+P2+P3;
fimplicit(f)

My expression looks like a sum of three parts P1, P2 and P3. fimplicit always draws a straight line like y=P_in. I see, that Mathlab writes that x and y values are "1x1 sym". I think, there is a mistake. My x should changes from 0 to 20 with some step (step doesn't matter). I tryed to write "x=0:0.1:20;" before the expression for f(x,y), but Mathlab wrote: "Invalid indexing or function definition. Indexing must follow MATLAB indexing. Function arguments must be symbolic variables, and function body must be sym expression". Also I tryed to write in expression for fimplicit: fimplicit(f [0 20 0 10]) but it only caused a changing of axes, but the graph still looked like a line.

How can I fix this problem? The graph should be a curve, not a straight line

Answer 1

[Asvin Kumar]

Try plotting with the fsurf command as follows:

fsurf(f, [-10 10 -10 10])
fsurf(f, [-10 10 4 6]) % or this 

You’ll notice that the z value of the function is 0 along the line y=5. This is also the reason why you see a straight line when using fimplicit. It plots the solution to f(x,y)=0 which in this case is y=5.

You can verify this by evaluating the following expressions:

vpa(subs(f,[x y],[4 5])) 
vpa(subs(f,[x y],[-1e4 5])) 
vpa(subs(f,[x y],[1e4 5]))

These three expressions evaluate the function value at (4,5) , (-10000, 5) and (10000, 5) respectively. The large values should give you an indication of how slowly the function value deviates from 0 along y=5.