Solving this system of differential equation with NDSolve

Question

I'm trying to solve this system of differential equation (the ones in yellow markers are just constants):

The initial conditions for this system are n(0) = NN, na(0) = Na, n'(0) = q(0) = 0, q'(0) = some_constants, na'(0) = Ip(0)×Some_other_constants (shown in my code below).

FWHM = 9;
sd = 2.355*FWHM;
\[Sigma]p = 1.15*10^-21; (* m^2 *)
\[Sigma]pa = 3.8*10^-21; (* m^2 *)
\[Sigma]e11 = 1.1*10^-20;  (* m^2 *)
\[Sigma]111 = 0.4*10^-20;  (* m^2 *)
\[Sigma]aa1 = 4.1*10^-20;  (* m^2 *)
\[Sigma]ea1 = 0.9*10^-20;  (* m^2 *)
\[Sigma]1a1 = 1.3*10^-20;  (* m^2 *)
\[Sigma]eaa = 3.05*10^-20;  (* m^2 *)
\[Sigma]1aa = 2.0*10^-20;  (* m^2 *)

NN = 9*10^27; (* m^-3 *)
Na = 1.08*10^26; (* m^-3 *)

c = 3*10^8; (* m/s*)
\[Eta] = 1.32;
\[CapitalOmega] = 5.2*10^-8;
\[CapitalOmega]a = 2.5*10^-6;
L = 0.003; (*m*)
V = 0.38;
\[Tau] = 4*10^-9 ;(*s*)
\[Tau]a = 4*10^-9; (*s*)

Ip [t] ==  2.5 E^(-(t - 8)^2/(2 sd))*10^28;
   
sol = NDSolve[{
    n'[t] == Ip[t] \[Sigma]p (NN - n[t]) - (\[Sigma]e11 c)/\[Eta] n[t] q[t] -
      n[t]/\[Tau], 
    q'[t] == ((\[Sigma]e11 - \[Sigma]111) c)/\[Eta] n[t] q[t] - 
      q[t]/\[Tau]c[t] + (\[CapitalOmega] n[t])/\[Tau] - 
      (\[Sigma]aa1 - c)/\[Eta] na[t]q[t] + (\[Sigma]ea1*c)/\[Eta] (Na - na[t]) q[t] - 
      (\[Sigma]aa1 c)/\[Eta] (Na - na[t]) q[t], 
    na'[t] == -Ip[t] \[Sigma]pa na[t] + (\[Sigma]ea1 c)/\[Eta] (Na - na[t]) q[t] + 
      (Na - na[t])/\[Tau]a - (\[Sigma]aa1 c)/\[Eta] na[t] q[t] + 
      (2 \[CapitalOmega]a)/(\[Tau]a L (\[Sigma]eaa - \[Sigma]1aa))*
       ((exp ((\[Sigma]eaa - \[Sigma]1aa) (Na - na[t]) L)) - 1)  , 
    \[Tau]c [t_] == (\[Eta] L^3)/(8 c \[Pi]^2) (n[t] (\[Sigma]e11 - \[Sigma]111) V)^2, 
    Ip [t_] ==  2.5 E^(-(t - 8)^2/(2 sd))*10^28, 

    n[0] == NN , q[0] == 0, na[0] == Na, 
    n'[0] == 0, q'[0] == (\[CapitalOmega] NN)/\[Tau], 
    na'[0] == -Ip[0] \[Sigma]pa Na, 
    \[Tau]c[0] == (\[Eta] L^3)/(8 c \[Pi]^2) (NN (\[Sigma]e11 - \[Sigma]111) V)^2},
   {n, q, na}, {t, 0, 17}];

Plot[Evaluate[{n[t], na[t], q[t], \[Tau][t]} /. s], {t, 0, 17}, 
 PlotPoints -> (t, 0, 1000)]

I've tried to write my code as shown above, however I got an error message

NDSolve::idelay: Initial history needs to be specified for all variables for delay-differential equations.

I'm not familiar with this topic and I personally think that I am not intending to solve a 'delay' equation. Can anyone pointing out my mistakes? Thank you.

Answer 1

I fixed the exp that @cvgmt noted, pulled τc and Ip out of NDSolve, and got rid of the extra initial conditions.

τc[t] := (η L^3)/(8 c π^2) (n[t] (σe11 - σ111) V)^2;
Ip[t] := 2.5 E^(-(t - 8)^2/(2 sd))*10^28;

eqns = {n'[t] == Ip[t] σp (NN - n[t]) - (σe11 c)/η n[t] q[t] - n[t]/τ, 
q'[t] == ((σe11 - σ111) c)/η n[t] q[t] - q[t]/τc[t] + (Ω n[t])/τ
- (σaa1 - c)/η na[t] q[t] + (σea1*c)/η (Na - na[t]) q[t]
- (σaa1 c)/η (Na - na[t]) q[t], 
na'[t] == -Ip[t] σpa na[t] + (σea1 c)/η (Na - na[t]) q[t] + (Na - na[t])/τa
- (σaa1 c)/η na[t] q[t] + (2 Ωa)/(τa L (σeaa - σ1aa))*(Exp[(σeaa - σ1aa) (Na - na[t]) L] - 1)};

sol = NDSolve[{eqns, n[0] == NN, q[0] == 0, na[0] == Na}, {n, q, na}, {t, 0, 17}];

which results in

NDSolve::nderr -- Error test failure at t == 1.740427851496955`*^-8; unable to continue.

To get a hint why that might happen, I evaluated the time-derivatives at t==0:

eqns /. t -> 0 /. {n[0] -> NN, q[0] -> 0, na[0] -> Na}
(* {n'[0] == -2.25*10^36, q'[0] == 1.17*10^29, na'[0] == -2.26701*10^33} *)

Those are some absolutely huge numbers, so it's not surprising that NDSolve can't deal with this over a time range of t==0 to t==17. Are you sure those parameter values are correct? (I guess if the variables are also of large magnitude, it might OK?)

ps. out of curiosity, what do the equations model?