i=1; for jk=0.00001:0.01:10 dk=@(deltaUk)((4)^(1/3)*bi^(1/3)*deltaUk^(2/3)*eps0^(1/3)*(y+1)^(1/3))/(jk^(1/3)); f2=@(deltaUk) (A_const.*p.*int(exp((B_const.*p.*(dk.^2))./(2.*deltaUk.*(x-dk))),x,0,dk)-log(1+(1./y))); Uk=fzero(f2,[0 3000]); Mass1(1,i)=jk; Mass2(2,i)=Uk; i=i+1; end disp(Mass1);