function dcoeff = ODEsystem_3F_DE(dt, coeff, x, uu, r, zeta)

% System of ODE that solves for the coefficients of the MGF/CF function
% INPUTS: dt     ---> time step (used by the matlab ODE solver)
%         coeffv ---> value of the coefficient
%         x      ---> vector with the model (2-factor) parameters
%         uu     ---> point where to compute the MGF(uu)
%         rf     ---> risk-free rate
%         d      ---> dividend yield
%
% OUTPUT: dcoeff ---> MGF coefficients (real and imaginary part):
%                     alpha_real
%                     alpha_img  
%                     beta1_real  
%                     beta1_img
%                     beta2_real
%                     beta2_img
%                      
% Nicola Fusari (2011):
% You can distribute and use the code. Please cite the author if you use
% the code in in your research.
                      
                      
% Parameters assignment
rho1   = x(1);
v1     = x(2);
k1     = x(3);
s1     = x(4);


rho2   = x(5);
v2     = x(6);
k2     = x(7);
s2     = x(8);

v3     = x(9);
k3     = x(10);
s3     = x(11);

c0m    = x(12);
c0p    = x(13);
c1m    = x(14);
c1p    = x(15);
c2m    = x(16);
c2p    = x(17);
c3m    = x(18);
c3p    = x(19);

lambdam  = x(20);
lambdap  = x(21);

muv1_m   = x(22);
muv1_p   = x(23);


muv3_m   = x(24);  
muv3_p   = x(25); % (not used)

rho_j    = x(26); 
rho_jp   = x(27); % (not used)

ss3 = x(28);

% Mean of the independed Jumps
ind_muv3_m = muv3_m*rho_j;

% Mean of the contemporaneous Jumps
dep_muv3_m = muv3_m - ind_muv3_m;


ls_flag = 1;

% Compute the correlated jumps compensator: (Double Check Here)
compm = DE_MGF(1,0,0, -1, lambdam, lambdap, muv1_m, muv1_p, dep_muv3_m, muv3_p, ls_flag) - 1;
compp = DE_MGF(1,0,0,  1, lambdam, lambdap, muv1_m, muv1_p, dep_muv3_m, muv3_p, ls_flag) - 1;


dcoeff = zeros(8,1);

ur = real(uu);
ui = imag(uu);


% Correlated jumps plus and minus
[c_phimr, c_phimi] = DE_MGF(uu, coeff(3)+1i*coeff(4), coeff(7)+1i*coeff(8),  -1, lambdam, lambdap, muv1_m, muv1_p, dep_muv3_m, muv3_p, ls_flag);
[c_phipr, c_phipi] = DE_MGF(uu, coeff(3)+1i*coeff(4), coeff(7)+1i*coeff(8),   1, lambdam, lambdap, muv1_m, muv1_p, dep_muv3_m, muv3_p, ls_flag);

c_phim = (c_phimr + 1i* c_phimi); 
c_phip = (c_phipr + 1i* c_phipi);

% Independent jumps plus and minus
[ind_phimr, ind_phimi] = DE_MGF(0.0, 0.0, coeff(7)+1i*coeff(8),  -1, lambdam, lambdap, 0.0, 0.0, ind_muv3_m, 0.0, ls_flag);

ind_phim = (ind_phimr + 1i* ind_phimi);

% Jump pos
phipr = real(c_phip);
phipi = imag(c_phip);

% Jump Neg
phimr  = real(c_phim+ind_phim -1);
phimi  = imag(c_phim+ind_phim -1);


% alpha

dcoeff(1) = r*ur + c0m*(-1 + phimr - compm*ur) + c0p*(-1 + phipr - compp*ur) + ...
            coeff(3)*k1*v1 + coeff(5)*k2*v2 + coeff(7)*k3*v3 - ur*zeta;

dcoeff(2) = c0p*phipi + c0m*(phimi - compm*ui) + coeff(4)*k1*v1 + coeff(6)*k2*v2 + coeff(8)*k3*v3 - ui*(c0p*compp - r + zeta);

% beta1

dcoeff(3) = (1/2)*(-2*coeff(3)*k1 - coeff(4)^2*s1^2 + coeff(3)^2*s1^2 - ui*(2*coeff(4)*rho1*s1 + ui) - ur +...
            ur*(2*coeff(3)*rho1*s1 + ur) + 2*c1m*(-1 + phimr - compm*ur) + 2*c1p*(-1 + phipr - compp*ur));

dcoeff(4) = c1m*phimi + c1p*phipi - (1/2)*(1 + 2*c1m*compm + 2*c1p*compp)*ui + ...
            ui*(coeff(3)*rho1*s1 + ur) + coeff(4)*(-k1 + s1*(coeff(3)*s1 + rho1*ur));


% beta2

dcoeff(5) = (1/2)*(-2*coeff(5)*k2 - coeff(6)^2*s2^2 + coeff(5)^2*s2^2 - ui*(2*coeff(6)*rho2*s2 + ui)...
            - ur + ur*(2*coeff(5)*rho2*s2 + ur) + 2*c2m*(-1 + phimr - compm*ur) + ...
            2*c2p*(-1 + phipr - compp*ur));
   
dcoeff(6) = c2m*phimi + c2p*phipi - (1/2)*(1 + 2*c2m*compm + 2*c2p*compp)*ui + ui*(coeff(5)*rho2*s2 + ur) ...
            + coeff(6)*(-k2 + s2*(coeff(5)*s2 + rho2*ur));


        
% beta3

dcoeff(7) = (1/2)*(-2*coeff(7)*k3 - coeff(8)^2*s3^2 + coeff(7)^2*s3^2 - 2*c3p*compp*ur + ...
            2*c3m*(-1 + phimr - compm*ur) - ss3*(ur + ss3*(ui - ur)*(ui + ur)));


dcoeff(8) = (-coeff(8))*k3 + c3m*phimi + coeff(8)*coeff(7)*s3^2 - c3m*compm*ui...
            - c3p*compp*ui - (ss3*ui)/2 + ss3^2*ui*ur;
        
        

        

