Hello all,
I am using solve to find the numerical value of two constants from two polynomial equations. As I vary the number of terms in the equation (set by changing the value of nmax), solve will return a solution for nmax <= 35, but not for higher values, instead stating: Warning: Could not extract individual solutions. Returning a MuPAD set object.. Do you have any suggestions on how I might be able to find values for my constants at higher polynomial orders?
Many thanks!
More info: I have analytically solved a differential equation, and the solution is the sum of two power series with two unknowns. For certain values of my constants, I can truncate the series after only a few terms and get a convergent solution. However, as I vary the value of my constants, I need to go to higher orders to get a convergent solution.
Code:
% Specify the appropriate constants for finding the bounds of the depletion region
k = 8.6173324E-5; % The Bolzmann constant (eV/K)
T = 300; % The temperature (K)
q = 1.60219E-19; % The fundamental electrical charge (C)
VT = k*T; % The thermal voltage (eV)
NA = 5E16; % The p-type, core doping (cm^-3)
ND = 9E17; % The n-type, shell doping (cm^-3)
NC = 5.758E17; % Effective density of states in the conduction band (cm^-3)
NV = 9.543E18; % Effective density of states in the valence band (cm^-3)
EG = 1.6071; % The material bandgap (eV)
ni = sqrt(NC*NV)*exp(-EG/(2*VT)); % The intrinsic carrier density (cm^-3)
Phi_0 = VT*log(NA*ND/ni^2); % The built in voltage (eV)
V = 0; % The applied voltage (V)
d1 = 0.1E-4; % Emitter thickness (cm)
d2 = 1.9E-4; % The core radius (cm)
eps_GaP = 11.1; % The dielectric constant of GaP
eps_GaAs = 12.9; % The dielectric constant of GaAs
eps0 = 8.854187817620E-14; % Vacuum permittivty in F/cm
eps = (11.1 + 1.8*0.9)*eps0; % The overall permittivity
% Solve for the limits of the depletion region
syms x2 x4;
S = solve(Phi_0 + V == q*NA/(6*eps)*d2^2+q/(3*eps*d2)*(NA*x4^3+ND*x2^3)+q /(2*eps)*(ND*x2^2-NA*x4^2), NA*(d2^3-x4^3) == ND*((d2+x2)^3-d2^3));
for a = 1:length(S.x2)
if isreal(S.x2(a)) && S.x2(a)>0 && isreal(S.x4(a)) && S.x4(a)>0
x2_val = double(S.x2(a));
x4_val = double(S.x4(a));
end
end
% Specify the appropriate constants for finding the currents in the quasi-neutral regions
nmax = 40; % Set the sum limit to a finite, maximum value
Ln = 1E-4; % The diffusion length in the n-type region
Lp = 1E-6; % The diffusion length in the p-type region
mun = 7661; % The electron mobility (cm^2/V/sec)
mup = 367.5; % The hole mobility (cm^2/V/sec)
Dn = mun*k*T; % The diffusion coefficient for electrons (cm^2/sec)
Dp = mup*k*T; % The diffusion coefficient for holes (cm^2/sec)
N_im = 0.1; % The imaginary part of the index of refraction
wl = 6E-5; % The wavelength, in cm.
alpha = 4*pi*N_im/wl; % The absorption coefficent (cm^-1)
inc_pow = 100; % AM1.5G power (mW/cm^2)
h = 6.626E-34; % Planck's constant (Jsec)
c = 3E10 ; % Speed of light (cm/s)
flux = inc_pow/1000*wl/(c*h); % The incident photon flux (photons/cm^2/sec)
% Solve for the current in the quasi-neutral shell
bvals = sym('bvals',[nmax 1]); % Define a matrix to hold the symbolic values for the polynomial coefficients
syms b0 B; % Define the constants to be solved for with the boundary conditions
% Populate the coefficient matrix
bvals(1,1) = b0;
bvals(2,1) = -alpha*b0;
bvals(3,1) = -1/6*(4*alpha*bvals(2,1)+(alpha^2-1/Lp^2)*bvals(1,1)+alpha*flux/Dp);
for i = 4:nmax
bvals(i,1) = 1/(i-1)/i*(2*alpha*(i-1)*bvals(i-1,1)+(alpha^2-1/Lp^2)*bvals(i-2,1));
end
% Calculate the various components of the minority carrier density
C2 = 0;
C3 = 0;
for k = 1:nmax
C2 = (d2+x2_val)^(k-1)*bvals(k,1) + C2;
C3 = (d1+d2)^(k-1)*bvals(k,1) + C3;
end
Hg2 = 0;
Hg3 = 0;
for l = 1:nmax
Hg2 = (1/(Lp^2*(l+1)*(l+2)))^(l-1)*(d2+x2_val)^(2*(l-1)) + Hg2;
Hg3 = (1/(Lp^2*(l+1)*(l+2)))^(l-1)*(d1+d2)^(2*(l-1)) + Hg3;
end
Sp = 100;
C2 = C2*exp(alpha*(x2_val-d1)) + B*Hg2;
C3 = C3 + B*Hg3;
dC3 = 0;
for p=1:nmax
dCa = (p-1)*bvals(p,1)*(d2+d1)^(p-2);
dCb = alpha*bvals(p,1)*(d2+d1)^(p-1);
dCc = B*2*(p-1)*(1/(Lp^2*(p+1)*(p+2)))^(p-1)*(d1+d2)^(2*p-3);
dC3 = dCa+dCb+dCc+dC3;
end
S3 = solve(C2 == ni^2/ND*(exp(V/VT)-1),-Dp*dC3==Sp*C3);
