This is my Matlab Script file. Whenever,I run the script, it shows runtime error. What is the error in this script? Thanks

close all
clear all
syms phi
E=80*1.6e-22;
hb=1.01e-34;
vf=1e6;
v0=200*1.6e-19;
m=9.1e-31;
e1=1.6e-19;
W=50e-9;
D=0:10:100;
kx=sqrt(E*(2*m/hb.^2)*cos(phi));
ky=sqrt(E*(2*m/hb.^2)*sin(phi));
qx=sqrt(((E-v0).^2./(hb.^2*vf.^2))-ky.^2);
%qx=sqrt(2*m*(E-v0)/hb.^2)*cos(phi);
thta=atan(ky./qx);
bta=sqrt((E-v0)*(2*m/hb.^2)*(1+sin(thta).^2));
z=((kx.^2+bta.^2).^2./(4*kx.^2*bta.^2));
T=(1+z*sinh(bta*D).^2);
Tb=1./T;
G0=e1^2*m*vf*W/hb^2;
%Tb=(16*E/v0)*exp(-2*qx*D);
G=G0*int(cos(phi)*Tb ,-pi/2, pi/2);
plot(D,G/G0)
title('Conductivity as a fxn of barrier width D for bilayer');
xlabel('barrier width');
ylabel('conductivity');

Answer 1

[Walter Roberson]

Your D is a vector. Your D is defined in terms of D so T is a vector. Your Tb is defined in terms of T so Tb is a vector. cos(phi)*Tb is therefore a vector. You are trying to int() a vector. That might be a problem, but possibly not. You would be safer doing the elements one by one.

More of a problem is that your expressions to be integrated are messy, and probably closed form integrals do not exist. You might want to switch to numeric integrals. If you do, you will need to do the elements one by one for sure.

Your integral requires high precision to evaluate at all meaningfully; you are not going to be able to evaluate it using integral(). When I evaluate at 5000 digits, I still see spikes as if there were singularities in the denominator -- but algebraically there are no singularities in the denominator. I am checking now about oddities in the numerator.

[Walter Roberson]

The expression involves sin() of an expression that includes sin(theta) as a term. The function turns out to have a high frequency component to it, and a number of large terms of varying sign in the order of 10^98 to 10^108. Catastrophic cancellation if not done at high precision. The sinh() component turns out to be mostly imaginary, but that is then squared, resulting in a negative number. At regular intervals the denominator becomes orders of magnitude less than the numerator, resulting in a narrow pulse. I have experimentally determined that the pulse is maximum where the imaginary component of the sinh() term goes through 0, indicating a real-valued solution of the sinh() at isolated points -- and the real component is 0 there. When 0 is substituted in for the sinh() term, there is a singularity in the integration. I have shown that the sinh() term does cross 0.

The conclusion has to be that there integral as a whole has many singularities and so cannot be integrated. The singularities are narrow so it may be difficult to find them unless you plot over a sufficiently narrow interval -- an interval on the order of 1E-12 wide shows multiple periods and you really need to get down to 1E-15 wide to figure out where the singularity is.

Answer 2

[Pritesh Shah]

In symbiolic, you can not plot until you simply that expression for different values.

[Roshan]

So, what should I do now ?