kill(all)$ /* 初期化 */
t_max:50$ /* 走査の時間 */
div:1000$ /* 分割数 */
fpprec:5$ /* 有効数字桁数 */
I_p0:1*10^24$ /* 励起光の強さの最大値 */
dt:t_max/div$
I(t):=I_p0*exp(-(t-tmax/2)^2/0.6^2)$ /* 寿命0.6nsの励起光パルス t=t_max/2で最大になるよう調整 */
array(N,div)$array(Q,div)$array(P,div)$
N[0]:0$Q[0]:0$P[0]:0$
for s:0 while s Ndot:bfloat(I(s*dt))*2.4b-17*(3b18-N[s])-3.020833333b-6*N[s]*Q[s]-N[s]/5b-9,
Qdot:(1.2b3*N[s]+3.020833333b-6*N[s]*Q[s])/2-Q[s]/2b-10,
N[s+1]:N[s]+bfloat(Ndot*dt*1.0b-9),
Q[s+1]:Q[s]+bfloat(Qdot*dt*1.0b-9),
P[s+1]:Q[s+1]*0.75b6
)$
data:makelist([s,P[s]],s,0,div)$
plot2d([discrete,data]);
数式処理ソフトMaxima(Windows版)を使って。
これを書くのに無駄に五日間もかかった。
疲れた。