Problem in drawing curves after solving ode equations

Hi there.

I've written the following program.This program is part of a laser program

main program

tspan = [0 3e-9];                   % time interval, up to 2 ns    
[t,y] = ode45(@rate_eq_program_1,tspan,y0);
size(t);
t=t*1e9;
param_rate_eq            % input of needed parameter
r_2=1-(abs(r_R)).^2;
plot(t,r_2)

function

function ydot = rate_eq_program_1(t,y)
%
 param_rate_eq              % input of needed parameters             
%
 current1 = 18d-2;                % bias current (step function) [A]; 400mA 
 

sigmaa=(2.*epsilon0.*ref_index.*ref_indexg)./(hbar.*w_s).*((1+(abs(r_R)).^2).*(1-r_R))./(conf.*V_g.*g_n.*(y(1)-N_0));

parameter

c = 3d10;                  % velocity of light [cm/s]
e = 1.6021892d-19;         % elementary charge [C]
h_Planck = 6.626176d-34;   % Planck constant [J s]    
hbar = h_Planck/(2.0*pi);  % Dirac constant [J s]
% Geometrical dimensions
L = 5d-4;             % cavity length [cm];                  5 um
w = 9d-5;               % active region width [cm];            900 nm
d = 80d-8;              % thickness of an active region [cm];  80 Angstroms 
%V = L*w*d;  
V_a = 5.26d-7;              % volume of active region[cm^3] --------paper 2016
conf = 0.5;            % confinement factor [dimensionless]
conf_NC =0.3;          % Nanocavity confinement factor [dimensionless]
V_m = V_a/conf;           % cavity volume [cm^3]
V_NC=2.4d-7;              %Nanocavity volume[cm^3]
ref_index = 3.5;        % effective mode index
ref_indexg=3.5;       % Group refractive index
V_g = c/ref_index;      % group velocity [cm/s]
tau1 = 5d-10;          % Waveguide carrier lifetime and Nanocavity carrier lifetime [s]
g_n = 5d-16;           % differential gain (linear model) [cm^2]
N_0=1d18;               % carrier density at transparency [cm^-3]
N_ss=0.3d18;               % carrier density at steady state [cm^-3]
%current_th = 184d-3;   % current at threshold [A]; 184 mA
henry_i=100000;         %Internal loss factor[cm^-1]  (alphai)
p=1;                     % cavity parity with P=1(?1)
epsilon0=8.85*10^-14;        % free space permittivity[F.cm^-1]
rho=(2*epsilon0*ref_index*c)./gamma_c.*hbar.*w_r;  %normalisation factor
r_L=1;                            % Left mirror re?ectivity
%delta_w=(w_c-w_s);
delta_w=0.52*gamma_T;

The problem is that r_R is a function of y(2) and it is also noteworthy that the other curves obtained are correct. the value r_R remains unchanged while y(2) has changes.I expect the curve plot(t,r_2) to be as follows.(black curve)

But to my surprise, I get the following curve.

What can be the problem?

Thanks in advance for your help.

Best regards.

Answer 1

[ME]

I think your issue comes from the

y(:,2)-N_0

part of your r_R expression. This is asking to take a number away from a number and I doubt MATLAB is using enough decimal places to be able to capture such a large subtraction from such a small number.

Unfortunately, right now I don't have a fix for you but I'm pretty certain that is your problem. Perhaps somebody else has a suggestion on how to address that issue?

[mohammad heydari]

I know this problem but I don't know how to fix it