MATLAB: Differential equations of time with two parameters given as functions of time.

MATLAB C/C++ Graphics Libraryode45

I have tried to solve the following differential equations of time. When I put Omega_C and Omega_F ,which are two parameters, as functions of time, I obtained errors. I look forward to your kind help..
function MTMA1(theta_M,theta_G,R_GM,gamma,time_final)
el_c=1.602e-19;hbar=(6.6261e-34)/(2*pi);T=300;gamma_M=1e14;
epsilon_0=8.854e-12; k_B=1.381e-23; E_F=1.36*el_c; c=3e8;
lambda=570e-9; omega=2*pi*c/lambda; omega_p=omega; omega_c=omega_p;
L_x=0.5e-9; d=L_x; L_z=7e-9; L_y=L_z; R_M=15e-9;
V_G=L_x*L_y*L_z; %V_G=pi*L_z^2*L_x;
V_M=(4/3)*pi*R_M^3;
phi_1=pi-0.5*pi-theta_G; phi_2=pi-0.5*pi-theta_M; theta_Q=phi_1+phi_2;
%-------------------------------





mu_12=1e-29;mu_13=mu_12;
gamma_12=1.36e-6*el_c/hbar; gamma_13=gamma_12; gamma_32=0.35*gamma_12;
%-------------------------------
epsilon_b=12.9; epsilon_q=6.5; epsilon_inf=5.7;
epsilon_star=(2*epsilon_b+epsilon_q)/(3*epsilon_b);
epsilon_G= 1 + ...
(1i*el_c^2/(8*epsilon_0*hbar*d*omega))*el_c* ...
(tanh((hbar*omega + 2*E_F)/(4*k_B*T)) + ...
tanh((hbar*omega - 2*E_F)/(4*k_B*T))) + ...
(el_c^2/(8*pi*hbar*epsilon_0*d*omega))*el_c* ...
log(((hbar*omega+2*E_F)^2)/((hbar*omega-2*E_F)^2+(2*k_B*T)^2))-...
(el_c^2/(pi*hbar*epsilon_0*d*omega))*el_c*(E_F/(hbar*omega + ...
1i*hbar*gamma));
epsilon_M=epsilon_inf-((omega_p)^2/((omega)^2+...
1i*omega*gamma_M));
%-------------------------------
zeta_x=1 - pi*L_x/(2*L_z); zeta_z=pi*L_x/(4*L_z); %zeta_y= zeta_z;
alpha_Gx=4*pi*V_G*(epsilon_G-epsilon_b)/(3*epsilon_b + ...
3*zeta_x*(epsilon_G-epsilon_b));
alpha_Gz=4*pi*V_G*(epsilon_G-epsilon_b)/(3*epsilon_b + ...
3*zeta_z*(epsilon_G-epsilon_b));
alpha_M=V_M*((epsilon_M-epsilon_b)/(epsilon_M+2*epsilon_b));
%-------------------------------
% For the distances in the system:
% Note that R_GM is a variable;
R_QG=(sin(theta_M)/sin(theta_Q)*R_GM);
R_QM=(sin(theta_G)/sin(theta_Q)*R_GM);
%--------------------------------
pi_x=(1/(4*pi*epsilon_star))*(alpha_Gx*(3*cos(phi_1-1))/R_QG^3 +...
alpha_M* (3*cos(phi_2-1))/R_QM^3);
pi_z=(1/(4*pi*epsilon_star))*(alpha_Gz*(3*cos(theta_G-1))/R_QG^3 +...
alpha_M* (3*cos(theta_M-1))/R_QM^3);
%-------------------------------
phi_x=(-1*alpha_Gx*alpha_M/((R_GM^3)*(4*pi*epsilon_star)^2))*...
((3*cos(phi_1)-1)/R_QG^3 + 3*cos(phi_2)-1/R_QM^3);
phi_z=(2*alpha_Gz*alpha_M/((R_GM^3)*(4*pi*epsilon_star)^2))*...
((3*cos(theta_G)-1)/R_QG^3 + 3*cos(theta_M)-1/R_QM^3);
%-------------------------------
Lambda_x=(mu_12^2/((4*pi*epsilon_star)^2*(hbar*epsilon_0*epsilon_b)))*...
(alpha_Gx*(3*cos(phi_1)-1)^2/R_QG^6 + alpha_M*(3*cos(phi_2)-1)^2/R_QM^6);
Lambda_z=(mu_13^2/((4*pi*epsilon_star)^2*(hbar*epsilon_0*epsilon_b)))*...
(alpha_Gz*(3*cos(theta_G)-1)^2/R_QG^6 + alpha_M*(3*cos(theta_M)-1)^2/R_QM^6);
omega_12=2.172436*(el_c/hbar);
omega_13=2.172432*(el_c/hbar);
Delta_p=omega_12-omega_p; Delta_c=omega_13-omega_c;
Delta_2=Delta_p-Delta_c;
%-------------------------------------------------------------------------%













% %



% The Fields parameters %
% %
%-------------------------------------------------------------------------%
% The interaction time b/w the QDs & the field:
t=0:1e-12:time_final;
D_FP=100e-15; D_FC=30e-15; t_p=1e-15; tau=30e-15; t_prime=30e-15;
Omega_c=(3/(sqrt(2*pi*D_FP^2)))*exp(-((t_prime-t_p+tau)^2)/(2*D_FP^2));
Omega_p=(3/(sqrt(2*pi*D_FC^2)))*exp(-((t_prime-t_p)^2)/(2*D_FC^2));
%-------------------------------------------------------------------------%
[t,rho]=ode45(@EqOfMotions,t,[0;0;0;0;1;0;0;0;0]);
plot(t./1e-12,imag((rho(:,2))),'bla');
function MEs=EqOfMotions(t,rho)
%-------------------------------------------------------------------------%
% %
% Equations of motion for the system %
% %
%-------------------------------------------------------------------------%
MEs(1,1)=-(gamma_12+gamma_13)*rho(1)+1i*(Omega_c*(pi_z+phi_z)+...
Lambda_z*rho(3))*rho(7)+...
1i*(Omega_p*(pi_x+phi_x)+...
Lambda_x*rho(2))*rho(4)-...
1i*(Omega_c*conj(pi_z+phi_z)+...
conj(Lambda_z)*rho(7))*rho(3)-...
1i*(Omega_p*conj(pi_x+phi_x)+...
conj(Lambda_x)*rho(4))*rho(2);
%-------------------------------------------------------------------------%
MEs(2,1)=-((0.5*gamma_13+0.5*gamma_12)+1i*(Delta_p-Lambda_x*(rho(5)-...
rho(1))))*rho(2)+...
1i*Omega_p*(pi_x+phi_x)*(rho(5)-...
rho(1))+...
1i*(Omega_c*(pi_z+phi_z)+...
Lambda_z*rho(3))*rho(8);
%-------------------------------------------------------------------------%
MEs(3,1)=-((0.5*gamma_13+0.5*gamma_12+0.5*gamma_32)+...
1i*(Delta_c-Lambda_z*(rho(9)-rho(1))))*rho(3)+...
1i*Omega_c*(pi_z+phi_z)*(rho(9)-...
rho(1))+...
1i*(Omega_p*(pi_x+phi_x)+...
Lambda_x*rho(2))*rho(6);
%-------------------------------------------------------------------------%
MEs(4,1)=-((0.5*gamma_13+0.5*gamma_12)-1i*(Delta_p-conj(Lambda_x)*(rho(5)-...
rho(1))))*rho(4)-...
1i*Omega_p*conj(pi_x+phi_x)*(rho(5)-...
rho(1))-...
1i*(Omega_c*conj(pi_z+phi_z)+...
conj(Lambda_z)*rho(7))*rho(6);
%-------------------------------------------------------------------------%
MEs(5,1)=gamma_12*rho(1)+gamma_32*rho(9)-1i*(Omega_p*(pi_x+phi_x)+...
Lambda_x*rho(2))*rho(4)+...
1i*(Omega_p*conj(pi_x+phi_x)+...
conj(Lambda_x)*rho(4))*rho(2);
%-------------------------------------------------------------------------%
MEs(6,1)=-(0.5*gamma_32-1i*Delta_2)*rho(6)-1i*(Omega_c*(pi_z+phi_z)+...
Lambda_z*rho(3))*rho(4)+...
1i*(Omega_p*conj(pi_x+phi_x)+...
conj(Lambda_x)*rho(4))*rho(3);
%-------------------------------------------------------------------------%
MEs(7,1)=-((0.5*gamma_13+0.5*gamma_12+0.5*gamma_12)-1i*(Delta_c-conj(Lambda_z)*(rho(9)-...
rho(1))))*rho(7)-...
1i*Omega_c*conj(pi_z+phi_z)*(rho(9)-...
rho(1))-...
1i*(Omega_p*conj(pi_x+phi_x)+...
conj(Lambda_x)*rho(4))*rho(8);
%-------------------------------------------------------------------------%
MEs(8,1)=-(0.5*gamma_32+1i*Delta_2)*rho(8)+1i*(Omega_c*conj(pi_z+phi_z)+...
conj(Lambda_z)*rho(7))*rho(2)-...
1i*(Omega_p*(pi_x+phi_x)+...
Lambda_x*rho(2))*rho(7);
%-------------------------------------------------------------------------%
MEs(9,1)=gamma_13*rho(1)-gamma_32*rho(9)-1i*(Omega_c*(pi_z+phi_z)+...
Lambda_z*rho(3))*rho(7)+...
+1i*(Omega_c*conj(pi_z+phi_z)+...
conj(Lambda_z)*rho(7))*rho(3);
%-------------------------------------------------------------------------%
end
end

Best Answer

First of all: remember about array/vector operators (vectorization)
Can't you just put Omega_C and Omega_F inside the function?