MATLAB: Numerical solutions of equation

numerical solution of error function equation

I'm trying to solve the following equation with respect to b
I need to find the value of b code is
clear all
E=1.5e8;
syms b
 a0=-4e7/400;
T0=600;
T=300;
temp=1/(2*a0*(T0-T));
X=5050700286489334129256038400000/(4349164374701853503175341475971*b) + exp(erfinv(erf(pi/b)/2 + (633825300114114700748351602688*((5879788993030477*E)/73786976294838206464 + 459358515080506015625/144115188075855872))/(4349164374701853503175341475971*b) - ((erf(pi/b) - (1267650600228229401496703205376*((5879788993030477*E)/73786976294838206464 + 459358515080506015625/144115188075855872))/(4349164374701853503175341475971*b))^2 + 1/5)^(1/2)/2)^2)*exp(-(2*b*erfinv(erf(pi/b)/2 + (633825300114114700748351602688*((5879788993030477*E)/73786976294838206464 + 459358515080506015625/144115188075855872))/(4349164374701853503175341475971*b) - ((erf(pi/b) - (1267650600228229401496703205376*((5879788993030477*E)/73786976294838206464 + 459358515080506015625/144115188075855872))/(4349164374701853503175341475971*b))^2 + 1/5)^(1/2)/2) - b^2*1i)^2/(4*b^2))*(50507002864893341292560384/(4349164374701853503175341475971*b) + (50507002864893341292560384*(erf(pi/b) - (1267650600228229401496703205376*((5879788993030477*E)/73786976294838206464 + 459358515080506015625/144115188075855872))/(4349164374701853503175341475971*b)))/(4349164374701853503175341475971*b*((erf(pi/b) - (101014005729786682585120768*E + 4040560229191467303404830720000000)/(4349164374701853503175341475971*b))^2 + 1/5)^(1/2))) + exp(erfinv(erf(pi/b)/2 + (633825300114114700748351602688*((5879788993030477*E)/73786976294838206464 + 459358515080506015625/144115188075855872))/(4349164374701853503175341475971*b) - ((erf(pi/b) - (1267650600228229401496703205376*((5879788993030477*E)/73786976294838206464 + 459358515080506015625/144115188075855872))/(4349164374701853503175341475971*b))^2 + 1/5)^(1/2)/2)^2)*exp(-(2*b*erfinv(erf(pi/b)/2 + (633825300114114700748351602688*((5879788993030477*E)/73786976294838206464 + 459358515080506015625/144115188075855872))/(4349164374701853503175341475971*b) - ((erf(pi/b) - (1267650600228229401496703205376*((5879788993030477*E)/73786976294838206464 + 459358515080506015625/144115188075855872))/(4349164374701853503175341475971*b))^2 + 1/5)^(1/2)/2) + b^2*1i)^2/(4*b^2))*(50507002864893341292560384/(4349164374701853503175341475971*b) + (50507002864893341292560384*(erf(pi/b) - (1267650600228229401496703205376*((5879788993030477*E)/73786976294838206464 + 459358515080506015625/144115188075855872))/(4349164374701853503175341475971*b)))/(4349164374701853503175341475971*b*((erf(pi/b) - (101014005729786682585120768*E + 4040560229191467303404830720000000)/(4349164374701853503175341475971*b))^2 + 1/5)^(1/2))) + (5050700286489334129256038400000*(erf(pi/b) - (1267650600228229401496703205376*((5879788993030477*E)/73786976294838206464 + 459358515080506015625/144115188075855872))/(4349164374701853503175341475971*b)))/(4349164374701853503175341475971*b*((erf(pi/b) - (101014005729786682585120768*E + 4040560229191467303404830720000000)/(4349164374701853503175341475971*b))^2 + 1/5)^(1/2)) + 7737125245533627/9671406556917033397649408
 vpasolve(X==temp,b)
but matlab(2016 version) gives the following error:
ans =
Empty sym: 0-by-1
equation also contain imaginary part which is written 1i in the expression of X
please help me for the same code for the calculation of X given below please have a look maybe this can help, code for the calculation of X is:
   clear all
syms  E b;
N0=1e6;
E0=1e8;
Eeff=-0.4e8;
a=19.3567;
p=0.1;
ed=0.4e-9;
efe=4.5*1e6*8.854e-12;
x=2*efe*(E-Eeff)/(a*b*p*sqrt(pi));
xm=erf(pi/b);
delta=0.05;
H_f=x-0.5*(x-xm+sqrt((x-xm).^2+4*delta));
theta=b*erfinv(H_f);
N1=N0*p;
J=N1*erf(pi/b);
p1=0.5*N1*exp(-b*b*0.5);
P_diel=2*ed*(E-E0);
p2=erf((2*theta+i*b*b)/(2*b))+erf((2*theta-i*b*b)/(2*b))+N1*erf(theta/b);%+P_diel;
P=p1+p2+P_diel;
X= diff(P,E);

Best Answer

%clear all
syms  E b;
N0=1e6;
E0=1e8;
Eeff=-0.4e8;
a=19.3567;
p=0.1;
ed=0.4e-9;
efe=4.5*1e6*8.854e-12;
x=2*efe*(E-Eeff)/(a*b*p*sqrt(pi));
xm=erf(pi/b);
delta=0.05;
H_f=x-0.5*(x-xm+sqrt((x-xm).^2+4*delta));
theta=b*erfinv(H_f);
N1=N0*p;
J=N1*erf(pi/b);
p1=0.5*N1*exp(-b*b*0.5);
P_diel=2*ed*(E-E0);
p2=erf((2*theta+1i*b*b)/(2*b))+erf((2*theta-1i*b*b)/(2*b))+N1*erf(theta/b);%+P_diel;
P=p1+p2+P_diel;
X= diff(P,E);
Eval=1.5e8;
a0=-4e7/400;
T0=600;
T=300;
temp=1/(2*a0*(T0-T));
eqn = subs(X, E, Eval) - temp;
B = linspace(-10,10);
Y = double( subs(eqn, b, B) );
plot(B, Y)
vpasolve(eqn, b, [6 7])
answer is approximately 6.688061088729162561767129937816
Related Question