How to model 5 coupled rate equations?

Hi, I'm trying to use dsolve to solve a system of 5 coupled rate equations. As a check I plot the total number of particles to make sure it remains constant throughtout. If phonon_dist is set to 0 this is what happens so the code works. But this removes one of the levels from system because it is only populated and de-populated by terms that include phonon_dist.

When phonon_dist is 0.5 all I see is N_ex decay, I'm seeing no other transient behaviour. I'm not sure why this is happening.

I have re-derived the rate equations multiple times and double checked the matrix is correct. From what I can see the terms all balance, negative terms are balanced by positive terms. So I'm not sure why the population of particles changes, I can't where I am losing them.

The solutions calculated by Matlab are very long. Is it that Matlab can't handle solutions this big? Or is it I've gone code blind and can't see an obvious error?

Thanks in advance.

if true
  % % Define lifetime variables
T_th = 10e-11;
T_cap = 10e-9; 
T_ph = 10e-20;
T_r2p = 10e-8;
T_nr = 10e-7;
T_rcb = 10e-8;

% define the number of excited phonons at temperature T phonon_dist = 5e-1;

% set up problems, represent u and v with symbolic functions syms N_ex(t) N_cb(t) N_2p(t) N_1s(t) N_gr(t)

% set up the differential equations. With N_ex, N_cb, N_2p, N_1s and N_gr

A = [ (-1/T_th) 0 0 0 0 ; (1/T_th) (-1/T_cap -2*phonon_dist/T_ph -1/T_rcb) (phonon_dist*1/T_ph) (phonon_dist*1/T_ph) 0 ; 0 (1/T_cap + phonon_dist*1/T_ph) (-2*phonon_dist/T_ph -1/T_r2p -1/T_nr) (phonon_dist*1/T_ph) 0 ; 0 (phonon_dist*1/T_ph) (phonon_dist*1/T_ph) (-2*phonon_dist/T_ph -1/T_nr) 0 ; 0 (1/T_rcb) (1/T_nr + 1/T_r2p) (1/T_nr) 0 ];

%()

B = [N_ex ; N_cb ; N_2p ; N_1s ; N_gr];

% set-up the differential eqn ordinary_diff_eqns = diff(B) == A*B;

% create the initial conditions Initial_conds = B(0) == [1000 ; 0 ; 0 ; 0 ; 0];

% solve the differential equations using above defined intial conditions [Soln_Struct] = dsolve(ordinary_diff_eqns, Initial_conds);

% extracts the solutions for individual levels. dsolve solves eqns then % puts the solutions in alphabetical order so solutions need extracted. N_ex_Soln = Soln_Struct(1).N_ex; N_cb_Soln = Soln_Struct(1).N_cb; N_2p_Soln = Soln_Struct(1).N_2p; N_1s_Soln = Soln_Struct(1).N_1s; N_gr_Soln = Soln_Struct(1).N_gr;

Total_electrons = N_ex_Soln + N_cb_Soln + N_2p_Soln + N_1s_Soln + N_gr_Soln;

% Plot the solutions to see the time dependence. %clf deletes from the current figure all graphics objects clf hold on grid on title('V6 - timescale = 5e-10') fplot(N_ex_Soln) xlim([0,5e-10]) ylim([0,1050]) xlabel('Time') ylabel('Population') fplot(N_cb_Soln) fplot(N_2p_Soln) fplot(N_1s_Soln) fplot(N_gr_Soln) fplot(Total_electrons) legend('N_{ex}','N_{cb}','N_{2p}','N_{1s}','N_{gr}','Total Electrons','Location','Best') end

[Karan Gill]

Try numeric solutions using a function like ode45. That might give you a better idea of what's going on.