Output is taking infinite loops.How to solve?

3 parts of the code are posted. Part-1 is Depletion, Part-2 is Depletion Region and Part-3 is Main. The functions depletion and depletion region are called in the main part of the code (Part-3).The error might be due to while loops in the depletion region(Part-2).

PART-1

function output = depletion(x3, V, x4, epsilonp, epsilonn,Na, Nd)

% depletion.m

%fzero(@depletion, x3old, [], Vbi - bias, 0, epsilonp, ...

% %is used to calculate the value of x3 that self - ...

%consistently gives the correct voltage drop across the ...

%pn junction, in the cylindrical geometry

 q = 1.60219e-19;
 d2 = x3 + x4
 x2 = sqrt(Na * (d2^2 - x4^2)/Nd + d2^2) - d2  % removed semicolon to check output value of x2;
 if x4 > 0
     %splitting expression for output into 5 terms to identify which term
     %results in complex value for output on 22/8/14 at 14:47
%      T1=V
%      T2=(q * Na/(4 * epsilonp) * (d2^2 - x4^2))
%      T3=q * Na/(2 * epsilonp) * x4^2 *log(d2/x4)
%      T4=q * Nd/(4 * epsilonn) * ((d2 + x2)^2 -d2^2)
%      T5= q * Nd/(2 * epsilonn) *(d2 + x2)^2 * (log((d2 + x2)/d2))
%      
 output = V - (q * Na/(4 * epsilonp) * (d2^2 - x4^2)- q * Na/(2 * epsilonp) * x4^2 *log(d2/x4) - q * Nd/(4 * epsilonn) * ((d2 + x2)^2 -d2^2) + q * Nd/(2 * epsilonn) *(d2 + x2)^2 * (log((d2 + x2)/d2)));

% removed semicolon to check output value of output;

%output = T1 - T2 - T3 - T4 + T5

elseif x4 == 0 output = V - (q * Na/(4 * epsilonp) * d2^2 - q *Nd/(4 * epsilonn) * ((d2 + x2)^2 - d2^2) +q * Nd/(2 * epsilonn) * (d2 + x2)^2 * (log((d2 +x2)/d2)))

% removed semicolon to check output value of output;

 end

%added this if statement to trap complex output on 24/8/14 at 11:28 am

 % if(any(imag(output)))
 % display('complex output');
 % display(output);
 % display(x2);
 % display(x3);
 % display(x4);
% output=0;

 end    
 %end

PART-2

function [x1, x2, x3, x4, Vbi] = depletionregion(Nd, Na, R,epsilonp, epsilonn, ni,bias, windowwidth)

% File name changed from A2.m to depletionegion.m by Ananth at 13/4/14 at

% 21:26 hrs

% changed x3old in the last argument to 1 on 22/8/14 14:36 pm% changed x3old in the last argument to 1 on

22/8/14 14:36 pm

% changed 1 back to x3old on 22/8/14 at 14:38 q = 1.60219E-19; kB = 1.3807e-23; T = 300; Vbi = ((kB * T) / q) * log(Na * Nd / ni^2)

if bias >= Vbi
x1 = 0;
x2 = 0;
x3 = 0;
x4 = 0;
Vbi = 0;
elseif bias<Vbi
x3guess = sqrt(2 * Nd * epsilonn * epsilonp * (Vbi- bias) / (q * Na * (Nd * epsilonn + Na * epsilonp)));
x3old = 0;
x3new = sqrt(2 * Nd * epsilonn * epsilonp * (Vbi -bias) / (q * Na * (Nd * epsilonn + Na * epsilonp)));
while abs(x3old - x3new)>1e-9
x3old = x3new;

x3new = fzero(@depletion, x3old, [], Vbi-bias, 0, epsilonp, epsilonn, Na, Nd);
    % changed x3old in the last argument to 1 on 22/8/14 14:36 pm
% changed 1 back to x3old on 22/8/14 at 14:38
  % trapping complex output from fzero function; 
  %if(x3new ==0)
      % break;
  % end

end
x3max = x3new;
x2max = sqrt(Na / Nd * x3max^2 + x3max^2) - x3max;

if (x2max + x3max) >= R
x4 = 0;
x1 = 0;
x2 = x2max * R / (x2max + x3max)
x3 = x3max * R / (x2max + x3max)
Vbi = (q * Na / (4 * epsilonp) * x3^2 - q * Nd / (4 * epsilonn) * ((x3 + x2)^2 - x3^2) +q * Nd / (2 * epsilonn) * (x3 + x2)^2 *(log((x3 + x2) / x3)));
elseif (x2max + x3max)<R
x4old = 1;
x4new = 0;

 while abs(x4old - x4new)>1e-9
 x4old = x4new

 x3 = fzero(@depletion, x3guess, [], Vbi -bias, x4old, epsilonp, epsilonn, Na, Nd);
 %trapping complex output from fzero function
 %if(x3new==0)
    % break;
 %end

 x2 = sqrt(Na * ((x3 + x4old)^2 - x4old^2) /Nd + (x3 + x4old)^2) - x3 - x4old;

 if x3 >= (R - windowwidth)
 x1 = R - x2 - x3;
 x4new = 0;
 elseif x3<(R - windowwidth)
x1 = max(0, windowwidth - x2);
 x4new = max(0, R - x1 - x2 - x3);
 end

 end
  x4 = x4new;

 else
 disp('error 2')
 end

 else
     disp('error 3')
  end

end

PART-3

%clears variables, closes windows and clears screen clear all; close all; clc;

%defining variables T= 310; % absolute temperature in Kelvin Dn=1.0; Dp=0.3; Nd=1e18; Na=1e18;

epsilonn=11.8; epsilonp=11.8; ni=1.5e10;

sigma_n=1e-5; sigma_p=1e-6;

% T1=0;T2=0;T3=0;T4=0;T5=0; Vth=0.026*T/300 ; window_width=1E-6;

alpha1 = 1.0;

%creating vectors len_start=0; len_end = 99; len_num_points =100;

vec_of_lengths = linspace(len_start,len_end, len_num_points);

rad_start=1; rad_end =100 ; rad_num_points = 100;

EField_num_points =100 ;

for count = 1:length(vec_of_lengths)

vec_of_rad = linspace(rad_start,rad_end, rad_num_points); vec__of_aspect_ratios = vec_of_lengths./vec_of_rad;

V_built_in=0.7;% built-in voltage in volts

volt_start = 0; volt_end = 0.95*V_built_in; volt_num_points = 100;

vec_of_volts = linspace(volt_start, volt_end, volt_num_points); length_vec_of_volts = length(vec_of_volts);

matrix1=[];

for count2= 1:17 matrix1(:,:,count2) = zeros(length_vec_of_volts, rad_num_points); end

for count2=1:4 matrix2(:,:,:,count2) = zeros(EField_num_points, rad_num_points, length_vec_of_volts); end

for count2=1:2 matrix3(:,count2) = zeros(1, rad_num_points); end

for count2=1:length(vec_of_rad)

for count3=1:length(vec_of_volts)

Ln=vec_of_rad(count3);% Ln is the diffusion length of electrons and is set equal to the radius.

Nr=Dn/(Ln*Ln*sigma_n*Vth);

taun0 = 1/(sigma_n*Nr*Vth);

taup0 = 1/(sigma_p*Nr*Vth);

taun=taun0;

taup=taup0;

trapdensitycount3(count3,count2)=taun0;

Lp=sqrt(taup*Dp);

bias(count3,count2) = vec_of_volts(count2); % assigning bias a value from voltage vector

%call depletion region function [x1,x2,x3,x4,V_built_in] = depletionregion(Nd, Na,vec_of_rad(count2),epsilonp, epsilonn, ni,bias(count3,count2), window_width)

if(bias(count3,count2) > V_built_in ) %V_built_in(count3,count2)) break; end %end if %add x1 and x2 d1=x1 + x2;

%add x3 and x4 d2=x3 + x4;

% define dimensionless variables

 beta1(count3,count2) = vec_of_rad(count2)/Lp;
 beta2(count3,count2) = (vec_of_rad(count2) - x1)/Lp;  % x1(count3,count2))/Lp; this commented part has to be used after 
 beta3 = alpha1*vec_of_lengths(count2); 
end

end

end

[Rick Rosson]

Please format your code so that it is readable. Make it easy for us to answer your question.

[John D'Errico]

This question is virtually impossible to answer by anyone who is not the author. We have no idea of the inputs to the code. No real idea of what it does or should do.

Instead, learn to use the debugger! Use it to trace through your code, watching what happens.

[Ced]

Try adding sprintf messages in your for loops to see where the problem is. Once you know that, add a line there saying "keyboard". This will stop the simulation at that place so you can have a look at all the values. Pressing F5 will then advance one iteration, so you can have a look at the same values etc.

[Guillaume]

while sprintf / keyboard / etc. would work, why not use the debugging tools designed in matlab. That is set breakpoints, step throught the code and stop on errors.