Plot a function with an integral involved

Dear all

I would like to plot the following function

where f(E-V,T) is

with μ set to zero , k_B*T = 0.0362 (unit : meV) and E replaced by E-V. So actually dependence on "T" is eliminated.

For the rest of the function

Here, delta is assumed to be 1.5 (meV) and Z = 1.

beta is as follow

I've tried so many methods but all of them have failed.

(Each part of the function is checked without errors)

when I drop the part of the derivativ of the function f, it works.

However, when I put it back, the result is trivial (without calculating the equation but returning somthing like int( ~~~~~, E, -3/2, 3/2) only)

Code

---------------------

syms E

V = [-5:0.1:5]

G = int ((exp((E-V)./ 0.362)./( 0.362*(exp((E-V)./ 0.362) + 1).^2)).*((2*(1+(E./(abs(1.5^2-E.^2)).^0.5).^2))./(((E./(abs(1.5^2-E.^2)).^0.5).^2)+(1+2*1^2)^2)), E, -1.5, 1.5) ...

+ int ((exp((E-V)./ 0.362)./( 0.362*(exp((E-V)./ 0.362) + 1).^2)).*((2*(E./(abs(1.5^2-E.^2)).^0.5))./(1+(E./(abs(1.5^2-E.^2)).^0.5)+2*1^2)), E, -inf, -1.5) ...

+ int ((exp((E-V)./ 0.362)./( 0.362*(exp((E-V)./ 0.362) + 1).^2)).*((2*(E./(abs(1.5^2-E.^2)).^0.5))./(1+(E./(abs(1.5^2-E.^2)).^0.5)+2*1^2)), E, 1.5, inf) ;

plot(V, G)

I did use integral instead of int to acquire the data but it only works when the variable V is dropped from the equation (Let V=zero)

Code

----------------------

fun = @(E)(exp((E)./ 0.362)./( 0.362*(exp((E)./ 0.362) + 1).^2)).*((2*(1+(E./(abs(1.5^2-E.^2)).^0.5).^2))./(((E./(abs(1.5^2-E.^2)).^0.5).^2)+(1+2*1^2)^2))

y = integral(fun, -1.5, 1.5)

Don't know how to sweep V while using the above code

Thanks for the help in advance