I am trying to plot a figure with different value of a parameter but the graph is plotting on negative Y axis.

I am trying to plot a figure with different value of 'a' but when I am entering value of 'a' below 0.5, the graph is plotting on negative axis. As the equation I am plotting is effective capapcity equation, it can't be negative. Same is happening when I am plotting 'fun'.

I have attached the figure for reference.

function y = difference(u,d1,n,a,m,T,PsByN_0,UmaxN_0)
d1=20;
n=10^-11.4;
ne=0.5;
m=2.7;
a=0.40;
T=1;
PsByN_0dB=15;                                                     
PsByN_0=10.^(PsByN_0dB/10);
UmaxdB = 15;
UmaxN_0=10.^(UmaxdB/10);
%u = 0.01678
u=0.001:0.001:0.5;

fun1 = @(u)  (1./u).*log(((d1^m)./(a.*ne.*PsByN_0.*T.*u+d1^m).*a)./(1-a));
fun2 = @(u) (1./u).*log(((-exp(u.*UmaxN_0).*(exp(-PsByN_0.*u)))./(u.*UmaxN_0+PsByN_0.*u)).*(PsByN_0.*u)-(PsByN_0.*u.*(exp(-PsByN_0.*u))).*(expint(u.*UmaxN_0+PsByN_0.*u))+(exp(-PsByN_0.*u))+((PsByN_0.*u).*(exp(-PsByN_0.*u))).*(expint(PsByN_0.*u))+(exp(u.*UmaxN_0))./((UmaxN_0/PsByN_0)+1))
fun = @(u) (fun1(u) - fun2(u));

%g0=fzero(fun,[0.001, 0.02])
g0 = fzero(fun,[0.01])

fun1(g0)
fun2(g0)

plot(u,fun1(u));
hold on;
grid on;

[dpb]

Something's wrong in your equations -- they never intersect for u>0 and there is a singularity at the origin so the soution is simply bogus.

>> [g0,fz,flg,o] = fzero(fun,[0.01]);
>> o.message
ans =
    'Current point x may be near a singular point. The interval [-0.0028, 0.019051] 
     reduced to the requested tolerance and the function changes sign in the interval,
     but f(x) increased in magnitude as the interval reduced.'
>> 

plotting the functions shows that fun1 approaches -inf as u-->0 while fun2 never crosses the origin.

You need to debug the equations and correct whatever is incorrect there or, if there really isn't a problem there, then change your expectations that there is an intersecting point.

[Walter Roberson]

log(-(447600772982339*a)/(137438953472*(5*10^(1/2)*a*u + 447600772982339/137438953472)*(a - 1)))/u - log(exp(-10*10^(1/2)*u) + exp(10*10^(1/2)*u)/2 + 10*10^(1/2)*u*exp(-10*10^(1/2)*u)*expint(10*10^(1/2)*u) - 10*10^(1/2)*u*exp(-10*10^(1/2)*u)*expint(20*10^(1/2)*u) - 1/2)/u

is your equation being plotted. For u=0.01 the a that gives 0 is about 0.504.

[dpb]

NB:

>> fun1
fun1 =
  function_handle with value:
    @(u)(1./u).*log(((d1^m)./(a.*ne.*PsByN_0.*T.*u+d1^m).*a)./(1-a))
>>

in which the argument of the log() function is

>> fnA
fnA =
  function_handle with value:
    @(u)(((d1^m)./(a.*ne.*PsByN_0.*T.*u+d1^m).*a)./(1-a))
>> 

Since log(x)<0 for x<1, any value of u for which fnA(u) < 1 is going to produce a negative value for fun1(u).

So,

>> fzero(@(x)fnA(x)-1,-10)
ans =
 -171.6444
>> 

shows you can't get to a positive value for fun1() until u is a large negative value. Meanwhile, that turns fun2() into complex plane.

[Dhawal Beohar]

So, is it my fun1 is incorrect or fun2?

[Walter Roberson]

Both might be incorrect, but at least fun2 is incorrect.

[Dhawal Beohar]

can I share my formulation with you? I am not sure where I am making the mistake...

[Walter Roberson]

You could post an image of the equations

[Dhawal Beohar]

I have attahced the formulkation for fun2. Please let me know if you find any mistake.

[dpb]

"Both might be incorrect, but at least fun2 is incorrect."

It's also possible both are correct but the hypothesis they intersect is wrong -- or one or more of the other constants is/are wrong to shift the solution for fnA() above to be >0 such that fun1() is also positive.

Or, they could be the correct solution to the wrong problem; or there are incorrect assumptions (or at least unfulfilled/inconsistent with) in the formulation of the problem itself.

We've no knowledge of the problem; no way to judge the applicability of the solution/methods and even if work through the math, certainly no idea of particular constants.