MATLAB: Problem solving laser rate equations with ode45 command

rate equations

Hi there, I have a question about modelling rate equations. Basically the modelling is solving differential equations hence I'm trying to use ODE45.The equations are as follows:
It should also be noted that the two values of and are as follows
I have written the above equations in a program as follows.
——Main program———–
clear all
tspan = [0 2d-9];                   % time interval, up to 2 ns
y0 = [0,0,0,0];     
[t,y] = ode45('rate_eq_program_1',tspan,y0);
size(t);
t=t*1d9;
y_max = max(y);
y1 = y_max(1);
y2 = y_max(2);
y3 = y_max(3);
y4 = y_max(4);
d = plot(t, y(:,1)/y1,'-.', t, y(:,2)/y2,'--', t, y(:,3)/y3,'-*', t, y(:,4)/y4,'-.');   % divided to normalize
xlabel('time [ns]','FontSize',14);  % size of x label
ylabel('Arbitrary units','FontSize',14);    % size of y label
set(gca,'FontSize',14);                 % size of tick marks on both axis
legend('carrier density','nanocavity carrier density','forward field','backward field')  % legend inside the plot
pause
close all
——-function——–
function y = rate_eq_program_1(t,y)
%

 param_rate_eq_fano              % input of needed parameters             
%
 current1 = 4d-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));     
ydot(1) = current1./(e.*V_a)-y(1)./tau1-(V_g.*g_n.*(y(1)-N_0).*sigmaa.*(abs(y(3))).^2)./V_m; 
ydot(2) = -y(2)./tau1-(conf_NC.*V_g.*g_n.*(y(2)-N_0).*rho.*(abs(y(4))).^2)./V_NC;
ydot(3) = 1/2.*(1-1i.*henry).*conf.*V_g.*g_n.*(y(1)-N_ss).*y(3)+gamma_L.*((y(4)./r_R-y(3)));
ydot(4) = (-1i.*delta_w-gamma_T).*y(4)-p.*gamma_c.*y(3)+1/2.*(1-1i.*henry).*conf_NC.*V_g.*g_n.*(y(2)-N_0).*y(4);
ydot = ydot';                   % must return column vector
——-param_rate_eq_laser——– 
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;          
w = 9d-5;            
d = 80d-8;            
%V = L*w*d;  
V_a = 5.26d-7;             
conf = 0.5;           
conf_NC =0.3;       
V_m = V_a/conf;         
V_NC=2.4d-7;              
ref_index = 3.5;       
ref_indexg=3.5;       
V_g = c/ref_index;      
tau1 = 5d-10;         
g_n = 5d-16;          
N_0=1d18;               
N_ss=3d18;               
henry_i=100000;       
henry=1;                
gamma_L=8.5*10^12;       
gamma_c=383*10^9;       
gamma_T=387*10^9;        
w_r=193*10^12;            
w_s=196*10^12;         
w_c=250*10^12;            
p=-1;                     
epsilon0=8.85*10^-14;        
rho=(2*epsilon0*ref_index*c)./gamma_c*hbar*w_r; 
r_L=1;                            
delta_w=(w_c-w_s);
And I get the following answer at the output.
0	0	0	0
0	0	0	0
0	0	0	0
0	0	0	0
0	0	0	0
0	0	0	0
0	0	0	0
0	0	0	0
0	0	0	0
0	0	0	0
0	0	0	0
0	0	0	0
0	0	0	0
0	0	0	0
0	0	0	0
0	0	0	0
0	0	0	0
0	0	0	0
0	0	0	0
0	0	0	0
0	0	0	0
0	0	0	0
0	0	0	0
0	0	0	0
0	0	0	0
0	0	0	0
0	0	0	0
0	0	0	0
0	0	0	0
0	0	0	0
0	0	0	0
0	0	0	0
0	0	0	0
0	0	0	0
0	0	0	0
0	0	0	0
0	0	0	0
0	0	0	0
0	0	0	0
0	0	0	0
0	0	0	0
Even if the given parameters have errors for any reason, this should not cause zero output.
Thanks in advance for helping!!

Best Answer

Change the function declaration line of your differential equation system function to: 
function ydot = rate_eq_fano_program_1(t,y)
That should work.
The problem is that it is returning ‘y’ (that is the input to the function), and is not returning ‘ydot’ that it should be returning.
Also, the problem with the plot call is that: 
y_max =
   232.9712e+012     0.0000e+000     0.0000e+000     0.0000e+000
so only ‘y(:,1)’ plots. The rest (all normalised by ‘ymax’) are Inf, and will not plot.
Related Question