I need to solve the following system of differential equations. enter image description here

When I have the solutions for $n_f$ and $v$, I need to find and plot $J=-e_\cdot n_{f} \cdot v$.

I wrote a code in matlab with all ODEs like this:

      function systemSolve

      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]; 


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

I suppose that:





I expect to find a current pulse like figure 1 but I get a exponential solution.

Figure 1

What is wrong in this code and can you suggest me how i can solve this problem?

  • $\begingroup$ Welcome to SciComp! Could you elaborate a bit on the physics of your model? Is "deltat" a problem parameter, or does that refer to a time step? $\endgroup$ Jun 29, 2013 at 20:35
  • $\begingroup$ Hi. deltat is a time step. And in the model the exponential factor for x'3 is e^((t/deltat)^2)) $\endgroup$
    – CarlosCr
    Jun 29, 2013 at 22:45
  • 1
    $\begingroup$ Having the time step for a numerical method inside your ODE system makes no sense; that information is for the discrete formulation only. $\endgroup$ Jun 30, 2013 at 19:48

1 Answer 1


You can use MATLAB's ode45, please refer to the first example here: Solve nonstiff differential equations; low order method.

Here I use $y$ instead of $x$. If you don't know how to write a $\mathtt{func.m}$ script and save it in your path, you can use function handle as well, in your case:

% output of function handle in vector form
rigid = @(t,y)[x(4); -x(2)+x(3);  ??? ; ???];
options = odeset('RelTol',1e-6,'AbsTol',[1e-6 1e-6 1e-6 1e-6]);
[t,y] = ode45(rigid,[0 T],[y10 y20 y30 y40],options);

??? is left for you to fill up, you can tweak the options if you read the document of the link above, y10 through y40 are the initial values. Lastly J = y(:,1).*y(:,3) then plot $J$ against $t$.

  • $\begingroup$ Hi. Thanks a lot for your assistance. I don´t have initial values in my problem. it is possible avoid this initial conditions?. or I have to deduce them. $\endgroup$
    – CarlosCr
    Jun 29, 2013 at 22:50
  • $\begingroup$ @user4671 You have to have initial condition, initial condition means what the system is like when you start your timer. $\endgroup$
    – Shuhao Cao
    Jun 29, 2013 at 23:21
  • $\begingroup$ Depending on the argument to the exponential, ode45 may be inappropriate due to stiffness, and replaced with ode15s. $\endgroup$ Jun 30, 2013 at 20:14
  • $\begingroup$ The matlab code that i wrote for this problem is the next: i´ve included some constants for simulation purpouses: $\endgroup$
    – CarlosCr
    Jun 30, 2013 at 21:03
  • 1
    $\begingroup$ en.m.wikipedia.org/wiki/Nondimensionalization $\endgroup$ Jul 1, 2013 at 2:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.