how to store the value of intersection points from fzero into a vector

I am calculating fzero using fzero function. For each value of 'a' (ranging from 0.1 to 1) I need to get value of u*(intersection of fun1 and fun2) as shown in below graph(just for reference). And then i want to store the value of u* in a different vector. Can someone please help me to get the values and store them.Below is my current code which is calculating u* only for a=0.4.

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.4;
T=1;
PsByN_0dB=5;                                                     
PsByN_0=10.^(PsByN_0dB/10);
UmaxdB = 5;
UmaxN_0=10.^(UmaxdB/10);
u=0.01:0.001:2;

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;
plot(u,fun2(u));
hold on;
grid on;

Answer 1

[Torsten]

Before using fzero or fsolve, you should check whether fun1 - fun2 really has a zero.

To do this, uncomment the lines

u=[0.01:0.01:2];

plot(u,fun(u,a(3)))

and insert the a-value (in this case a(3)) for which you want to perform the test.

For a>0.5, it seems that fsolve has problems to find a solution. Check it with the method I suggested.

d1=20;
n=10^-11.4;
ne=0.5;
m=2.7;
a=0.01:0.01:0.5;
T=1;
PsByN_0dB=5;                                                     
PsByN_0=10.^(PsByN_0dB/10);
UmaxdB = 5;
UmaxN_0=10.^(UmaxdB/10);
fun1 = @(u,a) (-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,a) (fun1(u,a) - fun2(u));
%u=[0.01:0.01:2];
%plot(u,fun(u,a(3)))
%g0=fzero(fun,[0.001,0.02])
options = optimset('Display','none');
g0 = arrayfun(@(a)fsolve(@(u)fun(u,a),[0.01],options),a)

g0 = 1×50

1.6739 1.4530 1.3231 1.2304 1.1581 1.0986 1.0480 1.0039 0.9646 0.9293 0.8971 0.8674 0.8399 0.8141 0.7900 0.7671 0.7455 0.7248 0.7050 0.6861 0.6678 0.6501 0.6330 0.6164 0.6001 0.5843 0.5688 0.5536 0.5386 0.5238

%fun(g0,a)
plot(a,g0)

[Dhawal Beohar]

Thank you for your time Torsten.

I am not getting that from where the value u is coming in for fun1 and fun2? when I am not uncommenting the lines you mentioned.

Also, what does these 50 values of g0 are? is it the u* which is calculaing i.e., the intersection of fun1 and fun2?

Please help.

[Torsten]

I am not getting that from where the value u is coming in for fun1 and fun2? when I am not uncommenting the lines you mentioned.

u is the variable in the function definition that "fsolve" tries to solve for - so it is a formal parameter that you don't need to give a value anywhere - except for the initial value 0.01 in the command

g0 = arrayfun(@(a)fsolve(@(u)fun(u,a),[0.01],options),a)

Also, what does these 50 values of g0 are? is it the u* which is calculaing i.e., the intersection of fun1 and fun2?

Yes, g0(i) is the u* corresponding to a(i) (i=1,...,50)

[Torsten]

@Dhawal Beohar comment moved here

Alright....So i don't need to use any for loop or something to iterate through the values of a and then store the value of u* in other vector? I think the solution you provided serves the purpose right. I am just being double sure, as this is really important for me right now. Thanks a ton !!!

[Torsten]

Alright....So i don't need to use any for loop or something to iterate through the values of a and then store the value of u* in other vector?

This is all done in the line

g0 = arrayfun(@(a)fsolve(@(u)fun(u,a),[0.01],options),a)

If you have more experience with for-loops: this code is equivalent to the one above, but you have more flexibility providing the initial guess for u*:

d1=20;
n=10^(-11.4);
ne=0.5;
m=2.7;
a=0.01:0.01:0.5;
T=1;
PsByN_0dB=5;                                                     
PsByN_0=10.^(PsByN_0dB/10);
UmaxdB = 5;
UmaxN_0=10.^(UmaxdB/10);
fun1 = @(u,a) (-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,a) (fun1(u,a) - fun2(u));
options = optimset('Display','none');
g00 = 0.01;
g0 = zeros(size(a));
for i=1:numel(a)
    g0(i) = fsolve(@(u)fun(u,a(i)),g00,options);
    g00 = g0(i);
end
plot(a,g0)

[Torsten]

@Dhawal Beohar comment moved here

Thanks for explaining things. I really appreciate your help. Thank you so much for being there !!!