Speeding up ODE for long time spans

Hey, I am trying to solve an no-stiff ODE, describing laser cavity with an external modulated injection. I would like to solve it for large time spans and hence, would like suggestion on optimizing my code as much as possible.

So as said before, I am solving a set of two coupled ODE, using ODE45. With an modulated input( which is a bessel filtered square wave). I wish to solve the differential equation for a span of about 10^6 times the time period of input function. I am attaching the code below. Using the profiler, I have tried to optimise as much as possible, by playing around with relative tolerances, absolute tolerances, clearing unnecessary variables etc.

 paramfile='laser_param';
eval(paramfile);
%global variable to be shared between this and the ODE function handle
global R1 A_fr del_om omega1 psi1 del_N kappa A_0 sq_fun t1 f;
%%Computing steady state parameters
A_fr= sqrt(( J-(gamma_n*N_th))/gamma_p);
R1=db2mag(-20);
psi1=-0.4;
kappa=k*R1;
r1=roots([1 -(2*(kappa/gamma_p)*A_fr*cos(psi1)) -A_fr^2 -(gamma_n*2*kappa*A_fr*cos(psi1)/(gamma_p*g)) ]);
A_0= r1(r1>=0.999*A_fr);
del_om=(-kappa*sqrt(1+alpha^2)*A_fr*sin(psi1+atan(alpha)))/(A_0*2*pi);
del_N=-2*kappa*A_fr/A_0/g*cos(psi1);%matrix elements
del_om=del_om*2*pi;
%%Setting ODE options
options =odeset('RelTol',1e-3,'AbsTol',[1e-2,1e-2]);
%%Frequency of input function 
omega1=1*1e9*2*pi;
%%Defining Time span, corresponds to 5 Cycles of the input signal
T=5*pi*2/omega1;
%%Creating a lookup table for the ODE solver
t1=linspace(0,T,5*35);
%bessel filtering
[b,a]=besself(2,15*omega1/2/pi);[bz,az]=impinvar(b,a,1/(t1(2)-t1(1)));
%the required square function
sq_fun=filter(bz,az,1.5*square(omega1*t1));%(m*square(omega1*t1)) tsmovavg((m*square(omega1*t1)),'e',15)
%interpolating to create a lookup table
f=@(xq)interp1(t1,sq_fun,xq);
%calling the ODE solver with steady state as initial conditions
[t,Z] = ode45(@leqn_solver_fun4,[0 T],[A_0*exp(1i*psi1) N_th+del_N],options);%

I have commented appropriately, please do comment/suggest where I could cut corners, so that when I am dealing with large time spans I can save considerable amount of computational time. I am also attaching the ode solver function and parameter file. So, far for the script written below, I have been able to optimize it to about 3.1 sec, for 5 input cycles.

[yassine chafiq]

Hi Arjun Lyer , i'm working on the optical in injection laser, i want to know if y resolve your problem? Thank's