getting an error could not find explicit solution and Vs1 = [ empty sym ]

%-------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*Vs/k1)))*Theta*nstar;
%syms Vs 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,Vs)
%Qsc=sqrt(2)*(Nb+Pb)*e*Ld*sqrt((cosh(ub+(e*Vs/k1)))/cosh(ub)-(e*Vs/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
Vs= solve(eqns, Vs1)
%pretty eqns

[Walter Roberson]

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