how to solve implicit equation as a function of Vs1

I am stuck in this since almost a week. Tried many things to solve this equation. could not get solution. Please help me and Walter sir please do not close the question!!! All helps are appreciated..

%-------Quantitative evolution of chemisorption Processes on %Metal oxide Semiconductor--------------

% ----------Discreption of the Programme-----
% This Matlab program calculate the relation between the total surface
% coverage area and the surface potential of the metal oxide
% semiconductors.The surface potential depends upon the gas pressure and
% temperature during the chemnisorption process.This program uses the
% Wolkienstien's theory of adsorption which is more adequate than the
% conventional Langmuir model.

% Variable used-----
% P- Pressure of the gas
% T-Temperature
% q0= Energy released during the weak adsorption
% nunegative= frequency of oscillation of the adions
% nunuetral= frequency of oscillation of the nuetral adsorbates
% M=  Molecular mass of the oxygen atom
% Eg= Energy band gap of the semiconductor
% Ga= degeneracy factor for the chemisorption induces states
% Ni= Intrinsic carrier concentration
% K= Boltazmann constant
% Vs= Surface potential
% Nd= Dopant concentration
% Nb= Bulk electron concentration
% Pb= Bulk hole concentration
% Ld= Debye length
% Qs= Surface charge density
% Qsc= Space charge density
%

p= 10^-6
T=300;
K=1.3807*10^-23;
Nd=10^24;
q0=0.1*1.6*10^-19;
M=32*1.67*10^-27;
Ga=2;
Nstar=10^19
es=5.4*8.85*10^-14;
e=1.6*10^-19;
x=0.21*1.6*10^-19  %(Ecb-Ef)=x
y=0.8*1.6*10^(-19)    %(Ecs-Eas)=y
Nc=2.5*10^25;
Nv=1.04*10^25;
Eg=2.42*1.6*10^-19;
M=32*1.67*10^-27;
k1=4.14*10^-21
k=1;
Vs=0.1
S0=10^-19;

M=32*1.67*10^-27;

nstar=10^25;
nunegative=10^13;
nunuetral=10^13
Ni=sqrt(Nc*Nv*exp(-Eg/(k1)))
Nb=Nd/2+sqrt(Nd^2/2+Ni^2)
Pb=Ni^2/Nb
ub=(k1)*log(Ni^2/Nb)
Ld=sqrt(es*(k1)/(e^2*Nd))
a=(1 + (1/Ga)*exp(y/k1)*exp(-(x/k1))*exp((e*Vs)/k1))
b=((k*S0)/(sqrt(2*pi*M*k1)*nunuetral))*exp(q0/k1)
c=(1 + (1/Ga)*(nunegative/nunuetral)*exp(-(x/k1))*exp((e*Vs)/k1))
Beta =(b/c)*a
Theta=(Beta*p)/(1+Beta*p)

Qs=-e/(1+(Ga*exp(-y/k1))*exp(x/k1)*exp(-(e*Vs1/k1)))*Theta*nstar;
%syms Vs1 Ga y k1 e x Theta nstar  
%eqns=-e/(1+(Ga*exp(-y/k1))*exp(x/k1)*exp(-(e*Vs/k1)))*Theta*nstar==0
%Vs= solve(eqns,Vs1)
%Qsc=sqrt(2)*(Nb+Pb)*e*Ld*sqrt((cosh(ub+(e*Vs1/k1)))/cosh(ub)-(e*Vs1/k1)*tanh(ub)-1);
%syms Vs1 Ga y k1 e x Theta nstar Nb Pb Ld ub
eqns=-e/(1+(Ga*exp(-y/k1))*exp(x/k1)*exp(-(e*Vs1/k1)))*Theta*nstar+sqrt(2)*(Nb+Pb)*e*Ld*sqrt((cosh(ub+(e*Vs1/k1)))/cosh(ub)-(e*Vs1/k1)*tanh(ub)-1)==0

Vs1= solve(eqns, Vs1)
Vs1= vpasolve(-e/(1+(Ga*exp(-y/k1))*exp(x/k1)*exp(-(e*Vs1/k1)))*Theta*nstar+sqrt(2)*(Nb+Pb)*e*Ld*sqrt((cosh(ub+(e*Vs1/k1)))/cosh(ub)-(e*Vs1/k1)*tanh(ub)-1)==0, Vs1,'random',true)

[Karan Gill]

Duplicate of https://www.mathworks.com/matlabcentral/answers/348913-explicit-solution-could-not-be-found

[Walter Roberson]

This does look like a duplicate of that question, but it has a different numeric solution, so something about the constants or the equations must be different.

Answer 1

[Walter Roberson]

You have

Vs1= solve(eqns, Vs1)

so after that statement, Vs1 is a symbolic numeric value. Then in your next statement,

Vs1= vpasolve(-e/(1+(Ga*exp(-y/k1))*exp(x/k1)*exp(-(e*Vs1/k1)))*Theta*nstar+sqrt(2)*(Nb+Pb)*e*Ld*sqrt((cosh(ub+(e*Vs1/k1)))/cosh(ub)-(e*Vs1/k1)*tanh(ub)-1)==0, Vs1,'random',true)

you use Vs1 just before 'random', so you are asking to solve a pair of numeric expressions, one of which is certainly non-zero (having just been obtained from solve(); vpasolve() applied to non-zero numeric expressions will fail.

When you re-use variable names in symbolic expressions it is much more common than not to run into these kinds of problems, where you get confused about the correspondence between symbolic name and value. Use different variable names to avoid the problem.

Note: the 'random' option requires R2014a or later.

Note: there are at least two different solutions, one positive and one negative. You can select which one you want by adding an assume()

assume(Vs1<0)

before you do the solve() or vpasolve()