Problem solving a differential equation with time dependent parameter.

I am trying to solve a system of rate equations, with time dependent ionization rate w(t) (my time dependent parameters in my case): my time vector is -10:0.1:10 so w(t) is also a vector and my equations are

dN0/dt=-W(t)*N0(t) and dN1/dt=W(t)*N0(t)

Think of it as a system with two levels, and the population at each level for time t is N0(t)and N1(t). An alternative approach for the equation would be to use the dN/dt=w(t)[Ng-N(t)] where Ng is a number, the initial population in ground state.

I am trying to solve this problem using the ode functions. However, using ode45 is not possible, since it takes a lot of time to compute and I cannot compute it and the other solvers don't give the "expected" solution.

I define my function as:

function dNdt=rate_equations(t,y,ft,f)

f = interp1(ft,f,t);

% dNdt=[-f.*y(1) ;   f.*y(1)];
dNdt=f.*[2.4*1e19-y];
end

the comment denotes another alternative approach and in the main script I calculate the solution as:

initial_cond=[0]; %initial_cond=[2.4*1e19;0]; when using the system pair
t_span=-10:0.1:10;
a=wt';
b=W_time';
[tt,NN]=ode23t(@(t,y)rate_equations(t, y, a, b), t_span, initial_cond);

My time dependent parameter vector W_time is defined as:

E_a=5.14*1e9 ;
omega_au=4.134*1e16;
eta0=376.6;
E_omega=sqrt(2*eta0*1e14); %Field of fundamental
E_omega2=sqrt(2*eta0*2e13);% field of second harmonic
theta=pi/2;
wt=-10:0.01:10 ;
W_time=4.*omega_au.*(E_a./abs(E_omega.*cos(wt)+E_omega2.*cos(2.*wt+theta))).*exp((-2/3).*(E_a./abs(E_omega.*cos(wt)+E_omega2.*cos(2.*wt+theta))))

I don't get the results I expected. In theory I should be getting something like this:

but with the plot saturating at some point. When I use the stiff solvers I get some very slight oscillations. The non stiff ones and ode45 are calculating for ever.

So my question is, what I am doing wrong? And how can I improve the computing time of the ode45 solver?

[Bilal sadiq]

Did you get your solution?