MATLAB: How to Solve Non Linear Electronutrility condition in Space region region in semiconductor? Using vpasolve it is showing [empty syms] error while a theortical solution exist.

Question

if true
  % code
end
r=40*10^-9;
T=500;
N=10^13;
K=1.3807*10^-23;
es=12*8.85*10^-12;
Nd=10^17;
Eg=3.6*1.6*10^-19;
e=1.6*10^-19;
x=0.21*1.6*10^-19;  %(Ecb-Ef)=x
y=1*1.6*10^-19;    %(Ecs-Eas)=y
%Eg=3.6*1.6*10^-19;
Nc=2.4154*10^24;
Nv=1.7959*10^25;
Ni=sqrt((Nc*Nv)*exp(-Eg/(K*T)))
Nb=Nc*exp(-(x/(K*T)))
Pb=Nv*exp(-(3.6*e/(K*T)))
ub=log(Nb/Ni)
Ld=sqrt(es*(K*T)/(e^2*Nd))
es=12*8.85*10^-12;
e=1.6*10^-19;
Ld=sqrt(es*(K*T)/(e^2*Nd));
%E=E0+1/6*(r/Ld)^2
%Qsc=sqrt(2)*(Nb+Pb)*e*Ld*sqrt(cosh(ub+E0+1/6*(r/Ld)^2/cosh(ub))-(E0+1/6*(r/Ld)^2)*tanh(ub)-1);
%syms r
%Nd*int(1-exp(E0+1/6*(r/Ld)^2),0,2.64*10^-22);
%E=[-50,50]
%F= Nd*int(1-exp*(E0+1/6*(r/Ld)^2),0,V)+ 4*pi*r^2*[Nt/1+2*exp(((y)/(K*T))+E0+1/6*(r/Ld)^2)]
%fsolve(@myfun,E)
%sqrt(2)*(Nb+Pb)*e*Ld*sqrt(cosh(ub+(E0+1/6*(r/Ld)^2/cosh(ub)))-(E0+1/6*(r/Ld)^2)*tanh(ub)-1)+ 4*pi*r^2*(N/1+2*exp(((-1.21*1.6*10^-19)/(K*T))+(E0+1/6*(r/Ld)^2)))=0
syms E0
vpasolve(sqrt(2)*(Nb+Pb)*e*Ld*sqrt(cosh(ub+(E0+1/6*(r/Ld)^2/cosh(ub)))-(E0+1/6*(r/Ld)^2)*tanh(ub)-1)+ 4*pi*r^2*(N/1+2*exp(((-1.21*1.6*10^-19)/(K*T))+(E0+1/6*(r/Ld)^2)))==0,E0)

Answer 1

There appears to be a minimum, but not a root.

To illustrate —

fcn = @(E0) (sqrt(2)*(Nb+Pb)*e*Ld*sqrt(cosh(ub+(E0+1/6*(r/Ld)^2/cosh(ub)))-(E0+1/6*(r/Ld)^2)*tanh(ub)-1)+ 4*pi*r^2*(N/1+2*exp(((-1.21*1.6*10^-19)/(K*T))+(E0+1/6*(r/Ld)^2))))
[E0s, fval] = fsolve(fcn, 1)
E0s =
          -34.966552734375
fval =
         0.624442049383109

With a complex initial estimate:

[E0s, fval] = fsolve(fcn, 1+1i)
E0s =
          -31.6877692741181 +      3.04633319691543i
fval =
          0.201062230775014 -    0.0783028987275935i