Array issue in numerical integration

Disclaimer: Brand New to Matlab; Apologies in advance for syntax or structure errors in my question.

I'm trying to solve the equation in the image. The integral needs to be solved numerically. Ultimately I'm trying to create a plot of Amin vs T for T between 1 and 1000. I keep getting dimensional errors, among other things, though. I'm not sure how to set it up so I can solve Amin for a range of T values to plot Amin vs T.

[For the image, note that I'm doing a 1D case so you can ignore the summation over i]

k=1.38e-23; %Boltzmann constant%
n=2.11e29; %Number density of atoms%
v=6134.8; %Speed of sound in YTaO4%
T=[1:1000]; %Temperature range in Kelvin%
hbar= 1.055e-34; %h-bar constant%
theta=v*(hbar/k)*(6*pi^2*n)^(1/3); %cutoff frequency for each polarization in Kelvin%
y = theta/T;
fun = @(x) (x^3*exp(x))/(e^x-1)^2;
I = integrate(fun,1,y);
Amin=(pi/6)^(1/3)*k*n^(2/3)*v*(T/theta)^2*I;
plot(T,Amin)

I know I need to get the sizes of T and Amin equal but don't know how to do this in Matlab. Thanks in advance for any help anyone can offer.

Answer 1

[Star Strider]

This runs:

k=1.38e-23; %Boltzmann constant%
n=2.11e29; %Number density of atoms%
v=6134.8; %Speed of sound in YTaO4%
T=[1:1000]; %Temperature range in Kelvin%
hbar= 1.055e-34; %h-bar constant%
theta=v*(hbar/k)*(6*pi^2*n)^(1/3); %cutoff frequency for each polarization in Kelvin%
y = theta./T;
fun = @(x) (x.^3.*exp(x))./(exp(x)-1).^2;
I = arrayfun(@(z)integral(fun,1E-4,z),y);
Amin=(pi/6)^(1/3)*k*n^(2/3)*v*(T./theta).^2.*I;
plot(T,Amin)

You may need to tweak it to get the result you want.

I substituted ‘1E-4’ for 0, but it still gives:

Warning: Infinite or Not-a-Number value encountered.

I leave that to you to sort.