function dcoeff = ODEsystem_2F_Gaussian(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 (2015):
% 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);

l0  = x(9);
l1  = x(10);
l2  = x(11);
myc1  = x(12);
syc1  = x(13);
mvc1  = x(14);
rhoc1 = x(15);


% Compute the correlated jumps compensator:
comp1 = (exp(myc1+0.5*syc1^2)/(1-rhoc1*mvc1)) - 1.0;


dcoeff = zeros(6,1);

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

d=zeta;



% alpha

dcoeff(1) = (-d)*ur + r*ur - l0*(1 + comp1*ur) + coeff(3)*k1*v1 + coeff(5)*k2*v2 + ...
            (exp(myc1*ur + (1/2)*syc1^2*(-ui^2 + ur^2))*l0*((-(-1 + coeff(3)*mvc1...
            + mvc1*rhoc1*ur))*cos(ui*(myc1 + syc1^2*ur)) - mvc1*(coeff(4) + rhoc1*ui)*sin(ui*(myc1 + syc1^2*ur))))/...
            (1 + mvc1*(coeff(3)^2*mvc1 + mvc1*(coeff(4) + rhoc1*ui)^2 - 2*rhoc1*ur + mvc1*rhoc1^2*ur^2 + 2*coeff(3)*(-1 + mvc1*rhoc1*ur)));

dcoeff(2) = (-(d + comp1*l0 - r))*ui + coeff(4)*k1*v1 + coeff(6)*k2*v2 + (exp(myc1*ur + (1/2)*syc1^2*(-ui^2 + ur^2))...
            *l0*(mvc1*(coeff(4) + rhoc1*ui)*cos(ui*(myc1 + syc1^2*ur)) - ...
            (-1 + coeff(3)*mvc1 + mvc1*rhoc1*ur)*sin(ui*(myc1 + syc1^2*ur))))/...
            (1 + mvc1*(coeff(3)^2*mvc1 + mvc1*(coeff(4) + rhoc1*ui)^2 - 2*rhoc1*ur + mvc1*rhoc1^2*ur^2 + 2*coeff(3)*(-1 + mvc1*rhoc1*ur)));
% beta1

dcoeff(3) = (1/2)*(-2*coeff(3)*k1 - 2*l1 - coeff(4)^2*s1^2 + coeff(3)^2*s1^2 - ...
            2*coeff(4)*rho1*s1*ui - ui^2 + (-1 - 2*comp1*l1 + 2*coeff(3)*rho1*s1)*ur + ur^2 +  ...
            (2*exp(myc1*ur + (1/2)*syc1^2*(-ui^2 + ur^2))* ...
            l1*((-(-1 + coeff(3)*mvc1 + mvc1*rhoc1*ur))* ...
            cos(ui*(myc1 + syc1^2*ur)) -  ...
            mvc1*(coeff(4) + rhoc1*ui)*sin(ui*(myc1 + syc1^2*ur))))/(1 +  ...
            mvc1*(coeff(3)^2*mvc1 + mvc1*(coeff(4) + rhoc1*ui)^2 - 2*rhoc1*ur + ... 
            mvc1*rhoc1^2*ur^2 +  ...
            2*coeff(3)*(-1 + mvc1*rhoc1*ur))));

dcoeff(4) = (1/2)*(ui*(-1 - 2*comp1*l1 + 2*coeff(3)*rho1*s1 + 2*ur) + ...
            2*coeff(4)*(-k1 + s1*(coeff(3)*s1 + rho1*ur)) +  ...
            (2*exp(myc1*ur + (1/2)*syc1^2*(-ui^2 + ur^2))* ...
            l1*(mvc1*(coeff(4) + rhoc1*ui)* ...
            cos(ui*(myc1 + syc1^2*ur)) - (-1 + coeff(3)*mvc1 +  ...
            mvc1*rhoc1*ur)* ...
            sin(ui*(myc1 + syc1^2*ur))))/(1 +  ...
            mvc1*(coeff(3)^2*mvc1 + mvc1*(coeff(4) + rhoc1*ui)^2 - 2*rhoc1*ur +  ...
            mvc1*rhoc1^2*ur^2 + 2*coeff(3)*(-1 + mvc1*rhoc1*ur))));
% beta2

dcoeff(5) = (1/2)*(-2*coeff(5)*k2 - 2*l2 - coeff(6)^2*s2^2 + coeff(5)^2*s2^2 - ...
            2*coeff(6)*rho2*s2*ui - ui^2 + (-1 - 2*comp1*l2 + 2*coeff(5)*rho2*s2)*ur + ur^2 + ...
            (2*exp(myc1*ur + (1/2)*syc1^2*(-ui^2 + ur^2))*l2*((-(-1 + coeff(3)*mvc1 + mvc1*rhoc1*ur))*...
            cos(ui*(myc1 + syc1^2*ur)) -  mvc1*(coeff(4) + rhoc1*ui)*sin(ui*(myc1 + syc1^2*ur))))/...
            (1 + mvc1*(coeff(3)^2*mvc1 + mvc1*(coeff(4) + rhoc1*ui)^2 - 2*rhoc1*ur + mvc1*rhoc1^2*ur^2 + ...
            2*coeff(3)*(-1 + mvc1*rhoc1*ur))));
   
dcoeff(6) = ui*(-(1/2) - comp1*l2 + coeff(5)*rho2*s2 + ur) + coeff(6)*(-k2 + s2*(coeff(5)*s2 + rho2*ur)) + ...
            (exp(myc1*ur + (1/2)*syc1^2*(-ui^2 + ur^2))*l2*(mvc1*(coeff(4) + rhoc1*ui)*cos(ui*(myc1 + syc1^2*ur)) ...
            - (-1 + coeff(3)*mvc1 + mvc1*rhoc1*ur)*sin(ui*(myc1 + syc1^2*ur))))/...
            (1 + mvc1*(coeff(3)^2*mvc1 + mvc1*(coeff(4) + rhoc1*ui)^2 - 2*rhoc1*ur + ...
            mvc1*rhoc1^2*ur^2 + 2*coeff(3)*(-1 + mvc1*rhoc1*ur)));
        
        
        
        

