algebric loop error in PV panel simulation... please help......

i have modeled the PV panel using a script file. the panel works when open circuit voltage and short circuit current are 0.665 and 5.75 respectively. although if i try to increase the open circuit voltage to say 42 and current to 7 , the circuit shows a algebraic loop error...the script file is given as an input to the current controlled source and the output voltage is fed as input.

please tell how to remove this algebraic loop error. function Ia=Iphoton(Va,G,TaC)

% Out: Ia = Module operating current (A) % In: Va = Module operating voltage (V) % G = insolation level(1G = 1000 W/m^2) % TaC = module operating voltage(celsius)

k = 1.381e-23; % boltzman constant q = 1.602e-19; % electron charge n = 1.2; % ideality factor Eg = 1.12; % band gap energy %Ns = 1000; % no. of series cell %Np=2; % no. of parallel cell TrK = 298; % reference temp. i.e. 25 degree celcius in kelvin Voc_TrK =42; % Voc : open circuit voltage at reference temp. across each series cell Isc_TrK = 7; % Isc: short circuit current at reference temp across each parallel cell. a = 3.5/1000; % current to temperature coefficient. TaK = 273 + TaC; % opearating current in kelvin %Vc = Va / Ns; % voltage acroos each cell in the array % calculation of Isc depending upon temp and a. %Isc = Isc_TrK * (1 + (a * (TaK - TrK))); % the photon current which is dependent upon the solar insolation level. %Iph = G * Isc; IL_T1 = Isc_TrK * G; IL = IL_T1 + a * (TaK-TrK); % all constants combined together Vt_TrK = n * k * TrK / q; % b = Eg * q/(n*k) b = Eg * q /(n * k); % calculation of diode current Ir_TrK = Isc_TrK / (exp(Voc_TrK / Vt_TrK) -1); Ir = Ir_TrK * (TaK / TrK)^(3/n) * exp(-b * (1 / TaK -1 / TrK)); % calculation of series resistence dVdI_Voc = -0.00985; % Take dV/dI @ Voc from I-V curve of datasheet Xv = Ir_TrK / Vt_TrK * exp(Voc_TrK / Vt_TrK); Rs = - dVdI_Voc - 1/Xv; % constants at arrays working temperature. Vt_Ta = n * k * TaK / q; % Ia = Iph - Ir * (exp((Vc + Ia * Rs) / Vt_Ta) -1) % f(Ia) = Iph - Ia - Ir * ( exp((Vc + Ia * Rs) / Vt_Ta) -1) = 0 % Solve for Ia by Newton's method: Ia2 = Ia1 - f(Ia1)/f'(Ia1) Ia=zeros(size(Va)); % Initialize Ia with zeros for j=1:5; Ia = Ia - (IL - Ia - Ir * ( exp((Va + Ia * Rs) / Vt_Ta) -1))... / (-1 - Ir * (Rs / Vt_Ta) * exp((Va + Ia * Rs) / Vt_Ta)); end end