y =0.145; p = 66; A = 3; B = 34; eps0 = 8.85e-12; bi = 5.18e-6; jk = 0.1; dk=@(deltaUk)(4)^(1/3)*bi^(1/3)*deltaUk.^(2/3)*eps0.^(1/3)*(y+1).^(1/3)/(jk^(1/3)); f = @(u, deltaUk) exp(B*p*dk(deltaUk)./(2* deltaUk*(u-1))); J = @(deltaUk) integral( @(u) f(u, deltaUk), 0,1, 'ArrayValued',true); f2 = @(deltaUk) p*A*dk(deltaUk).*J(deltaUk) - log(1+1/y); deltaUk = 10000:100:100000; plot(deltaUk, f2(deltaUk)); xlabel('\Delta U_k'); ylabel('f2'); Uk=fzero(f2,[1 10000000])