MATLAB: Solve a system of three ODE

ode45system

I have to solve this system.
diff(Np,t)==(Gamma*a0*vg*(n-n1)/(1+Ec*Np)-1/Tphot)*Np + Beta*Gamma*n/Tcarr;
diff(n,t)==I/(q*Vact)-n/Tcarr-vg*a0*(n-n1)/(1+Ec*Np)*Np ;
diff(F,t)==Alpha/2*(Gamma*a0*vg*(n-n1)-1/Tphot);
I know the value of these constants:
q = 1.6*10^-19;
Vact = 2.6110^-17;
ng = 4.2;
c = 3*10^8;
vg = c/ng ;
Tcarr = 1.46*10^-9 ;
a0 = 9.01*10^-20;
n1 = 1.36*10^24;
Gamma = 0.0561;
L = 250*10^-6;
nr = 3.6;
Beta = 6.46*10^-5;
Alpha = 5;
Alphai = 2300;
R = ((nr-1)/(nr+1))^2;
Alpham = (1/L)*log(1/R);
Alphauk = Alphai + Alpham;
gth = Alphauk/Gamma;
Tphot = 1/(Gamma*vg*gth);
Ec = 0;
I = 100*10^-3;
Conditions:
n(0) = 1.36*10^24;
Np(0)=0;
F(0)=0;
This system has to be solved for period of time t = [0 100]*10^-9.
I don't know, can I get function for n(t), Np(t) and F(t) in analytic form, but I need to plot these three function separately.

Best Answer

The Symbolic Math Toolbox cannot find an analytic solution for your system. The next best option is to create an anonymous function and integrate it numerically with ode15s:
syms Np(t) n(t) F(t) t T Y
q = 1.6E-19;
Vact = 2.61E-17;
ng = 4.2;
c = 3E+8;
vg = c/ng ;
Tcarr = 1.46E-9 ;
a0 = 9.01E-20;
n1 = 1.36E+24;
Gamma = 0.0561;
L = 250E-6;
nr = 3.6;
Beta = 6.46E-5;
Alpha = 5;
Alphai = 2300;
R = ((nr-1)/(nr+1))^2;
Alpham = (1/L)*log(1/R);
Alphauk = Alphai + Alpham;
gth = Alphauk/Gamma;
Tphot = 1/(Gamma*vg*gth);
Ec = 0;
I = 100*10^-3;
DE_sys = [diff(Np,t)==(Gamma*a0*vg*(n-n1)/(1+Ec*Np)-1/Tphot)*Np + Beta*Gamma*n/Tcarr;
diff(n,t)==I/(q*Vact)-n/Tcarr-vg*a0*(n-n1)/(1+Ec*Np)*Np ;
diff(F,t)==Alpha/2*(Gamma*a0*vg*(n-n1)-1/Tphot)];
[DE_sys_vf, Subs] = odeToVectorField(DE_sys);
DE_sys_fcn = matlabFunction(DE_sys_vf, 'Vars',{t Y});
DE_sys_fcn = @(t,Y)[-(Y(1).*6.435714285714287e-12-8.75257142857143e12).*Y(2)-Y(1).*6.849315068493151e8+2.394636015325671e34;(Y(1).*3.610435714285715e-13-9.81330604838636e11).*Y(2)+Y(1).*2.482232876712329e3;Y(1).*9.026089285714286e-13-2.45332651209659e12];
To understand what ‘Y(1)’ and the rest represent, see the ‘Subs’ variable, so Y(1)=Subs(1), and so for the rest.