我需要求解以下微分方程组。

当我有解决方案时$n_f$和$v$，我需要找到并绘制$J=-e_\cdot n_{f} \cdot v$.

我在 matlab 中编写了一个包含所有 ODE 的代码，如下所示：

      function systemSolve
      clc


      tr=50e-12;        % Recombination lifetime
      n0=1e11;         % Density of free carriers (Recordar que es 1e17 cm-3)
      tc=2e-12;         % Trapping time
      ts=30e-15;        % Carrier scattering time
      m=0.067*9.11e-31; % effective mass GaAs
      ev=8.854e-12;     % permitivity  
      n=900;            % factor geometry
      q=1.6e-19         % electron charge

      de=50e-15         %  delta t


      timeRange=[0 0.1e-12]; 
      initialConditionVector=[0;1e-15;1e-15;1e-15];
      [t,x]=ode45(@xprime,timeRange,initialConditionVector);
      figure(1),plot(t,x(:,1))

      J=q*x(:,3).*x(:,1);
      figure(2),plot(t,J)

      function f=xprime(t,x)
      f=[x(4); ...
        -(x(2)/tr)+(x(3)*q*x(1)); ...
        -(x(3)/tc)+(n0*(exp((t/de)^2))); ...
        -((1/ts)*x(4))-(((x(3)*q^2)/m*ev)*x(1)/n)+(q*x(2)/(m*n*ev*tr))];
      end
end

我想：

$x_1=v$

$x_2=P_{sc}$

$x_3=n_f$

$x_4=x'_1$

我希望找到一个如图 1 所示的电流脉冲，但我得到了一个指数解决方案。

这段代码有什么问题，你能建议我如何解决这个问题吗？