Code to numerically integrate a system of first-order ODEs

I need to solve the following system of differential equations.

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


I suppose that:

$x_1=v$

$x_2=P_{sc}$

$x_3=n_f$

$x_4=x'_1$

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

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

• 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? – Geoff Oxberry Jun 29 '13 at 20:35
• Hi. deltat is a time step. And in the model the exponential factor for x'3 is e^((t/deltat)^2)) – CarlosCr Jun 29 '13 at 22:45
• Having the time step for a numerical method inside your ODE system makes no sense; that information is for the discrete formulation only. – Geoff Oxberry Jun 30 '13 at 19:48

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$.

• 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. – CarlosCr Jun 29 '13 at 22:50
• @user4671 You have to have initial condition, initial condition means what the system is like when you start your timer. – Shuhao Cao Jun 29 '13 at 23:21
• Depending on the argument to the exponential, ode45 may be inappropriate due to stiffness, and replaced with ode15s. – Geoff Oxberry Jun 30 '13 at 20:14
• The matlab code that i wrote for this problem is the next: i´ve included some constants for simulation purpouses: – CarlosCr Jun 30 '13 at 21:03
• en.m.wikipedia.org/wiki/Nondimensionalization – Geoff Oxberry Jul 1 '13 at 2:25