v_B = sqrt((e*T_e)/(m_i*amu));
v_e = sqrt(8*k_B*T_ar/(pi*m_e));
L_pl = (l_B*m_e)/(e^2*n*A_sp);
A = (1/(2*e*epsilon*n*A_sp^2))
B = (1/(2*e*epsilon*n*A_sg^2))
f = @(t, x) [- x(4) - Jip + Jep*exp(- (e*A/T_e)*x(1)^2);
x(4) - Jig + Jeg*exp(- (e*B/T_e)*x(2)^2);
A*x(1)^2 - B*x(2)^2 + V0*cos(t) + x(3) - v_eff/Lpl_inv*x(4)];
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-8);
[t, xa] = ode45(f, [0:0.0001:0.06], [0, 0, -62, 0], options);
subplot(2,2,1), plot(t, xa(:,1)), grid on, xlabel('t'), title('x_{1}')
subplot(2,2,2), plot(t, xa(:,2)), grid on, xlabel('t'), title('x_{2}')
subplot(2,2,3), plot(t, xa(:,3)), grid on, xlabel('t'), title('x_{3}')
subplot(2,2,4), plot(t, xa(:,4)), grid on, xlabel('t'), title('x_{4}')
f = @(t, x) [- Jip + Jep*exp(- (e*A/T_e)*x(1)^2);
- Jig + Jeg*exp(- (e*B/T_e)*x(2)^2)];
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-8);
[t, xa] = ode45(f, [0:0.0001:0.1], [0, 0], options);
subplot(2,1,1), plot(t, xa(:,1)), grid on, xlabel('t'), title('x_{1}')
subplot(2,1,2), plot(t, xa(:,2)), grid on, xlabel('t'), title('x_{2}')