How to create array while using integral function?

Following is the code, there are parameters - b,c,gama1,sigma. And r that is the radius of the sphere is varing, also the time is varying. I want to get at each radial point at given time the value of I(integral) What am I doing is wrong?

K1=6;%(W/m.K)Thermal conductivity of Fe3O4 nanoparticle
K2=0.598; %Thermal conductivity of water-https://thermtest.se/thermal-conductivity-of-water
k1= 1.82E-06; %(m^2/s) Diffusivity of Fe3O4 nanoparticle
k2=2.29e-9; %at 25C (m^2/s) diffusivity of water- https://en.wikipedia.org/wiki/Molecular_diffusion
a=15e-9; %radius of nanoparticle
Vol= (4/3)*3.14*a^3;%Volume of nanoparticle
A=8.68E-05; %(W)heatproduced/ unit volume
b=(K2/K1)*((k1/k2)^1/2);
r=linspace(15e-9,17e-9,10);
c=1-(K2/K1);
t=linspace(1e-12,100e-12,10);
gama1=(a^2)/k1;
sigma=((r/a)-1)*((k1/k2)^1/2);

for t = 1e-12:100e-12
    r = 15e-9:17e-9;
    fun = @(y,gama1,sigma,b,c,t) ((exp(((-y.^2).*t)./gama1)).*(sin(y)-(y.*cos(y))).*((b.*y.*sin(y).*cos(sigma.*y))-((c.*sin(y)-y.*cos(y)).*sin(sigma.*y))))./((y.^3).*((((c.*sin(y))-(y.*cos(y))).^2)+((b.^2).*(y.^2).*(sin(y).^2))));
    I=integral(@(y) fun(y,gama1,sigma,b,c,t),0,inf); %integration wrt y with other constants
end

Answer 1

[VBBV]

K1=6;%(W/m.K)Thermal conductivity of Fe3O4 nanoparticle
K2=0.598; %Thermal conductivity of water-https://thermtest.se/thermal-conductivity-of-water
k1= 1.82E-06; %(m^2/s) Diffusivity of Fe3O4 nanoparticle
k2=2.29e-9; %at 25C (m^2/s) diffusivity of water- https://en.wikipedia.org/wiki/Molecular_diffusion
a=15e-9; %radius of nanoparticle
Vol= (4/3)*3.14*a^3;%Volume of nanoparticle
A=8.68E-05; %(W)heatproduced/ unit volume
b=(K2/K1)*((k1/k2)^1/2);
r=linspace(15e-9,17e-9,10);
c=1-(K2/K1);
t=linspace(1e-12,100e-12,10);
gama1=(a^2)/k1;
sigma=((r/a)-1)*((k1/k2)^1/2);
  
syms y

for J = 1:length(sigma)
    for K = 1:length(t)
        fun = ((exp(((-y.^2).*t(K))./gama1)).*(sin(y)-(y.*cos(y))).*((b.*y.*sin(y).*cos(sigma(J).*y))-((c.*sin(y)-y.*cos(y)).*sin(sigma(J).*y))))./((y.^3).*((((c.*sin(y))-(y.*cos(y))).^2)+((b.^2).*(y.^2).*(sin(y).^2))));
        I(K,J)= vpaintegral(fun,y,[0,Inf]); %integration wrt y with other constants
    end
end
I = (vpa(I,15))

I = 

Answer 2

[Dyuman Joshi]

%Formatting your code properly is a good coding practice

K1=6;                   %(W/m.K)Thermal conductivity of Fe3O4 nanoparticle
K2=0.598;               %Thermal conductivity of water-https://thermtest.se/thermal-conductivity-of-water
k1=1.82E-06;            %(m^2/s) Diffusivity of Fe3O4 nanoparticle
k2=2.29e-9;             %at 25C (m^2/s) diffusivity of water- https://en.wikipedia.org/wiki/Molecular_diffusion
a=15e-9;                %radius of nanoparticle
Vol=(4/3)*3.14*a^3;     %Volume of nanoparticle
A=8.68E-05;             %(W)heatproduced/ unit volume
b=(K2/K1)*((k1/k2)^1/2);
r=linspace(15e-9,17e-9,10);
c=1-(K2/K1);
t=linspace(1e-12,100e-12,10);
gama1=(a^2)/k1;
sigma=((r/a)-1)*((k1/k2)^1/2);

%preallocation
I = zeros(numel(t),numel(r));
%Function to integrate
fun = @(y,t) ((exp(((-y.^2).*t)./gama1)).*(sin(y)-(y.*cos(y))).*((b.*y.*sin(y).*cos(sigma.*y))-((c.*sin(y)-y.*cos(y)).*sin(sigma.*y))))./((y.^3).*((((c.*sin(y))-(y.*cos(y))).^2)+((b.^2).*(y.^2).*(sin(y).^2))));

for k = 1:numel(t)
    %Pass the value of t(k) 
    I(k,:)=integral(@(y) fun(y,t(k)),0,inf,'ArrayValued',true); %integration wrt y with other constants
end
format long
I

I = 10×10

5.253204730589755 5.253496876204714 5.253496836492238 5.253496764345236 5.253496718036804 5.253496699699072 5.253496651746186 5.253496649534816 5.253496645563740 5.253496581295285 5.250606466445789 5.253496875856392 5.253496824734424 5.253496798055983 5.253496781958984 5.253496757830739 5.253496678522883 5.253496630996006 5.253496573923628 5.253496627903624 5.248595986052435 5.253496881251547 5.253496844173969 5.253496807104954 5.253496770028961 5.253496733075693 5.253496695949829 5.253496658502814 5.253496621075591 5.253496585141805 5.246968241707677 5.253496883582442 5.253496838807810 5.253496794002702 5.253496781377496 5.253496746044602 5.253496684327441 5.253496648018087 5.253496617617959 5.253496636420552 5.245609714553339 5.253496881329098 5.253496844030694 5.253496806559629 5.253496770556284 5.253496733477224 5.253496695275669 5.253496658419976 5.253496621098609 5.253496589641397 5.244440063530290 5.253496881255291 5.253496844177262 5.253496807084363 5.253496770061161 5.253496732984319 5.253496695855780 5.253496658813329 5.253496621616801 5.253496585266910 5.243405230465747 5.253496881237234 5.253496844179250 5.253496807104066 5.253496770097128 5.253496736876893 5.253496707832046 5.253496604517689 5.253496711479302 5.253496482981266 5.242469707085853 5.253496881055251 5.253496844180170 5.253496807106646 5.253496770015861 5.253496733667430 5.253496699741413 5.253496643310109 5.253496648326808 5.253496553328025 5.241609991854334 5.253496879830237 5.253496844180530 5.253496807107732 5.253496770024883 5.253496733066472 5.253496697085797 5.253496654382595 5.253496629636052 5.253496574988831 5.240810143920073 5.253496874524455 5.253496844297906 5.253496807278024 5.253496770244370 5.253496733202896 5.253496696472889 5.253496657718357 5.253496624161022 5.253496581597618