Exponential Increase instead of Linear

The plot between V and I should be as shown, but it produces an exponential increase of I with V

q = 1.6e-19;
E = [-0.5:0.001:0.5];       % [eV]
V = [0:0.01:0.5];                % [V]
mu = 0;                           % [eV]
kB = 1.38e-23;                    % [J/K]
kB = kB/q;                        % [eV/K]
h = 6.626e-34/(44/7);             % [J.s]
h = h/q;                          % [eV.s]
T = 300;                          % [K]
g1 = 0.2;                         % [eV], gamma1
g2 = 0.2;                         % [eV], gamma2
g = 0.1;                          % [eV], g1*g2/g1+g2
dE = 0.001;                       % [eV], delta E

D = zeros(1,length(E));
mu1 = zeros(1,length(V));
mu2 = zeros(1,length(V));
f1 = zeros(1,length(V));
f2 = zeros(1,length(V));
f = zeros(1,length(V));
F = zeros(1,length(V));
I = zeros(1,length(V));
 
for j = 1:length(V)
     
   for i = 1:length(E)
    if E(i)>= -0.1
        D(i)=0.5;
    else
        D(i)=0;
    end
   end
 
    mu1(j) = mu + ((q*V(j))/2);
    mu1(j) = mu1(j)/q;                          % [eV]
    mu2(j) = mu - ((q*V(j))/2);
    mu2(j) = mu2(j)/q;                          % [eV]
    f1(j) = 1 ./ (1+exp((E(i)-mu1(j))/(kB*T)));
    f2(j) = 1 ./ (1+exp((E(i)-mu2(j))/(kB*T)));
    f(j) = f1(j)-f2(j);
    F(j) = D(i)*g*f(j)*dE;
    I(j) = ((2*q)/h)*sum(F(j));
    
end
 
plot(V,I,'b-','linewidth',2); hold on;
set(gca,'fontsize',20);
xlabel('V');
ylabel('I');
xlim([0 0.5]);
set(gca,'xtick',[0:0.1:0.5]);

[Looky]

Check your inner for loop. It computes the D values for index i 1:length(E) but only the i = length(E) is used afterwards. Also this loop does not contain the index j and E is nether changed, so why is it nested anyway? You could compute it beforehand?

Look at f1(j) = 1 ./ (1+exp((E(i)-mu1(j))/(kB*T))); Here the E(i) will always be the same value ( E(length(E))? is this intended?