Your E will not be a scalar, so having something "/" a calculation involving E would be invoking the "matrix right divide" operator, not the element-by-element division operator.
However even if you fix that, you have the issue that your function involves multiplying by all of q where q is a vector, so the integral would have to result in a vector (one for each element of q). Because integral() passes in a non-scalar for the parameter, that is incompatible -- unless you use the 'vectorized', true option which would pass in a scalar for E and so allow a vector result the same size as q to be returned. But... your context expects a scalar result.
Each time you include all of q in your fun you are using a q that has at least some of the entries still as their zeros() . Each of those leads to a zero result because the 0 will be multiplying the entire result of the term. But.. that includes the first time, where q is all zero, so the first output would return all 0, and no matter which of those all-zero values you selected to store into q(i), you would be storing a zero, so q would not change. And that implies that all of your results are going to come out zero.
I think you need to reconsider using q inside fun . Perhaps it was intended to be Q ?
Fs= 2.16*10^-5*pi; % Geometrical Factor for sun
Fa=pi; % Geometrical Factor for earth
Q=1.6*10^-16; % Electronic Charge
h= 4.136*10^-15; % Plancks Constant
c= 3*10^8; % Speed of light
K = 8.61*10^-5; % Boltzmanns Constant
Ts=5760; % Temparature of the sun
Ta=300; % Ambient Temparature
Eg=1.6; % BandGap Energy
A=((2*Fs)/(h^3*c^2));
B=(K*Ts);
C=((2*Fa)/(h^3*c^2));
D=(K*Ta);
vdata=0:0.01:1.6;
n=length(vdata);
q=zeros(1,n);
for i=1:n
V=vdata(i);
fun=@(E)((q.*E.^2).*((A./(exp(E./B)-1))-(C./(exp(E-(Q.*V))./B)-1)-(C./(exp(E./D)-1))));
whos
q(i)=integral(fun,Eg,inf);
end
plot(V,q)
