Vectorization a for loop which consists a transcendental function

Complete code:

tic
%Constants
Imp = 7.61;
Vmp = 26.3;
Pmax_e = 200.143;
Isc_n = 8.21;
Voc_n = 32.9;
Kv = -0.1230;
Ki = 0.0032;
Ns = 54;
k = 1.3806503e-23;
q = 1.60217646e-19;
T_kelvin = 273.15;
T_n = 25;            %STC Nominal Sicakligi 25C
G_n = 1000;          %STC Nominal Radiation 1000 W/m2

%Options
G = 1000;  %Gercek ortamin radyasyonu
T = 25;    %Gercek ortamin sicakligi
a = 1.3;
Vt = (Ns*k*(T+T_kelvin))/q;

%Block 1
dT = T - T_n;
Io = (Isc_n + Ki*dT)/(exp((Voc_n + Kv*dT)/(a*Vt))-1); % Eq. (7)
Rs = 0;
Rp_min = (Vmp/(Isc_n-Imp))- ((Voc_n-Vmp)/Imp); % Eq. (11)
Rp = Rp_min;


%Block 2
ep_max=1;
tol = 8e-5;
dV = 0.001;
while(ep_max>tol)
	Ipv_n = ((Rp+Rs)/Rp)*Isc_n;  % Eq. (10)
	Ipv = (Ipv_n + Ki*dT)*(G/G_n); % Eq. (4)
	Rp = (Vmp*(Vmp+Imp*Rs))/(Vmp*Ipv-Vmp*Io*exp(((Vmp+Imp*Rs)*q)/(Ns*a*k*(T+T_kelvin)))+Vmp*Io-Pmax_e); % Eq. (9)
	
	I_eski=Isc_n; %(V=0,I=Isc) point
	i=1;
	for V=0:dV:Voc_n
		I = I_eski;
		J = 1+((Io*Rs)/(Vt*a))*exp((V+Rs*I)/(Vt*a))+(Rs/Rp);
		F = I-Ipv+Io*(exp((V+Rs*I)/(Vt*a))-1)+((V+Rs*I)/Rp); % Eq. (3)
		I_yeni = I_eski - pinv(J)*F;
		P(i) = I_yeni*V;
		i=i+1;
		I_eski = I_yeni;
	end
	
	Pmax = max(P);
	ep_max = abs(Pmax-Pmax_e)
	Rs = Rs+0.0001;
end

fprintf("Ipv = %f\nIo = %g\na = %f\nRs = %f\nRp = %f\n", Ipv, Io, a, Rs, Rp);
toc

When I choose dV = 0.001 it takes too much time to complete. I've seen using Vectorization is very usefull for getting rid of for loops and it significantly reduces the execution time. But in my example the next values of I vector is determined by solving I_yeni = I_eski - pinv(J)*F equation. I don't have a complete values of I vector in the first place. Therefore how I can vectorize this for loop? Or is there an other fast way to calculate this values?

I have tried this:

dV = 0.001;
V_vector = 0:dV:Voc_n;
I_vector = zeros(1, length(V_vector)+1);
P = zeros(1, length(V_vector)); % pre-allocoating
while(ep_max>tol)
	Ipv_n = ((Rp+Rs)/Rp)*Isc_n;  % Eq. (10)
	Ipv = (Ipv_n + Ki*dT)*(G/G_n); % Eq. (4)
	Rp = (Vmp*(Vmp+Imp*Rs))/(Vmp*Ipv-Vmp*Io*exp(((Vmp+Imp*Rs)*q)/(Ns*a*k*(T+T_kelvin)))+Vmp*Io-Pmax_e); % Eq. (9)
	
	I_vector(1) = Isc_n; %(V=0,I=Isc) point
    for n=1:length(V_vector)
		J = 1+((Io*Rs)/(Vt*a))*exp((V_vector(n)+Rs*I_vector(n))/(Vt*a))+(Rs/Rp);
		F = I_vector(n)-Ipv+Io*(exp((V_vector(n)+Rs*I_vector(n))/(Vt*a))-1)+((V_vector(n)+Rs*I_vector(n))/Rp); % Eq. (3)
		
        I_vector(n+1) = I_vector(n) - pinv(J)*F;
        P(n) = I_vector(n+1)*V_vector(n);
    end
	Pmax = max(P);
	ep_max = abs(Pmax-Pmax_e)
	Rs = Rs+0.0001;
end

but in the end I still need to use a for loop to determine the I_vector's next values.