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)
%fun(g0,a)
plot(a,g0)
