## Problem in drawing curves after solving ode equations

on 11 Nov 2019 at 8:22
Latest activity Edited by mohammad heydari

on 12 Nov 2019 at 7:11

### ME (view profile)

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?
Best regards.

### ME (view profile)

on 11 Nov 2019 at 9:18
Edited by ME

### ME (view profile)

on 11 Nov 2019 at 9:23

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?