function [dcoeff1,dcoeff2,dcoeff3,dcoeff4,dcoeff5,dcoeff6] = ODEsystem_2F_Gaussian_gpu(h, dt, coeff1, coeff2, coeff3, coeff4, coeff5, coeff6, rho1, v1, k1, s1, rho2, v2, k2, s2, c0, c1, c2, muy, sy, mv, rhoj, 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
l0  = c0;
l1  = c1;
l2  = c2;
myc1  = muy;
syc1  = sy;
mvc1  = mv;
rhoc1 = rhoj;


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


dcoeff1=0;dcoeff2=0;dcoeff3=0;dcoeff4=0;dcoeff5=0;dcoeff6=0;

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

d=zeta;



% alpha

dcoeff1 = (-d)*ur + r*ur - l0*(1 + comp1*ur) + coeff3*k1*v1 + coeff5*k2*v2 + ...
            (exp(myc1*ur + (1/2)*syc1^2*(-ui^2 + ur^2))*l0*((-(-1 + coeff3*mvc1...
            + mvc1*rhoc1*ur))*cos(ui*(myc1 + syc1^2*ur)) - mvc1*(coeff4 + rhoc1*ui)*sin(ui*(myc1 + syc1^2*ur))))/...
            (1 + mvc1*(coeff3^2*mvc1 + mvc1*(coeff4 + rhoc1*ui)^2 - 2*rhoc1*ur + mvc1*rhoc1^2*ur^2 + 2*coeff3*(-1 + mvc1*rhoc1*ur)));
dcoeff1 = h*dcoeff1;

dcoeff2 = (-(d + comp1*l0 - r))*ui + coeff4*k1*v1 + coeff6*k2*v2 + (exp(myc1*ur + (1/2)*syc1^2*(-ui^2 + ur^2))...
            *l0*(mvc1*(coeff4 + rhoc1*ui)*cos(ui*(myc1 + syc1^2*ur)) - ...
            (-1 + coeff3*mvc1 + mvc1*rhoc1*ur)*sin(ui*(myc1 + syc1^2*ur))))/...
            (1 + mvc1*(coeff3^2*mvc1 + mvc1*(coeff4 + rhoc1*ui)^2 - 2*rhoc1*ur + mvc1*rhoc1^2*ur^2 + 2*coeff3*(-1 + mvc1*rhoc1*ur)));
dcoeff2 = h*dcoeff2;

% beta1

dcoeff3 = (1/2)*(-2*coeff3*k1 - 2*l1 - coeff4^2*s1^2 + coeff3^2*s1^2 - ...
            2*coeff4*rho1*s1*ui - ui^2 + (-1 - 2*comp1*l1 + 2*coeff3*rho1*s1)*ur + ur^2 +  ...
            (2*exp(myc1*ur + (1/2)*syc1^2*(-ui^2 + ur^2))* ...
            l1*((-(-1 + coeff3*mvc1 + mvc1*rhoc1*ur))* ...
            cos(ui*(myc1 + syc1^2*ur)) -  ...
            mvc1*(coeff4 + rhoc1*ui)*sin(ui*(myc1 + syc1^2*ur))))/(1 +  ...
            mvc1*(coeff3^2*mvc1 + mvc1*(coeff4 + rhoc1*ui)^2 - 2*rhoc1*ur + ... 
            mvc1*rhoc1^2*ur^2 +  ...
            2*coeff3*(-1 + mvc1*rhoc1*ur))));
dcoeff3 = h*dcoeff3;

dcoeff4 = (1/2)*(ui*(-1 - 2*comp1*l1 + 2*coeff3*rho1*s1 + 2*ur) + ...
            2*coeff4*(-k1 + s1*(coeff3*s1 + rho1*ur)) +  ...
            (2*exp(myc1*ur + (1/2)*syc1^2*(-ui^2 + ur^2))* ...
            l1*(mvc1*(coeff4 + rhoc1*ui)* ...
            cos(ui*(myc1 + syc1^2*ur)) - (-1 + coeff3*mvc1 +  ...
            mvc1*rhoc1*ur)* ...
            sin(ui*(myc1 + syc1^2*ur))))/(1 +  ...
            mvc1*(coeff3^2*mvc1 + mvc1*(coeff4 + rhoc1*ui)^2 - 2*rhoc1*ur +  ...
            mvc1*rhoc1^2*ur^2 + 2*coeff3*(-1 + mvc1*rhoc1*ur))));
dcoeff4 = h*dcoeff4;

% beta2

dcoeff5 = (1/2)*(-2*coeff5*k2 - 2*l2 - coeff6^2*s2^2 + coeff5^2*s2^2 - ...
            2*coeff6*rho2*s2*ui - ui^2 + (-1 - 2*comp1*l2 + 2*coeff5*rho2*s2)*ur + ur^2 + ...
            (2*exp(myc1*ur + (1/2)*syc1^2*(-ui^2 + ur^2))*l2*((-(-1 + coeff3*mvc1 + mvc1*rhoc1*ur))*...
            cos(ui*(myc1 + syc1^2*ur)) -  mvc1*(coeff4 + rhoc1*ui)*sin(ui*(myc1 + syc1^2*ur))))/...
            (1 + mvc1*(coeff3^2*mvc1 + mvc1*(coeff4 + rhoc1*ui)^2 - 2*rhoc1*ur + mvc1*rhoc1^2*ur^2 + ...
            2*coeff3*(-1 + mvc1*rhoc1*ur))));
dcoeff5 = h*dcoeff5;
   
dcoeff6 = ui*(-(1/2) - comp1*l2 + coeff5*rho2*s2 + ur) + coeff6*(-k2 + s2*(coeff5*s2 + rho2*ur)) + ...
            (exp(myc1*ur + (1/2)*syc1^2*(-ui^2 + ur^2))*l2*(mvc1*(coeff4 + rhoc1*ui)*cos(ui*(myc1 + syc1^2*ur)) ...
            - (-1 + coeff3*mvc1 + mvc1*rhoc1*ur)*sin(ui*(myc1 + syc1^2*ur))))/...
            (1 + mvc1*(coeff3^2*mvc1 + mvc1*(coeff4 + rhoc1*ui)^2 - 2*rhoc1*ur + ...
            mvc1*rhoc1^2*ur^2 + 2*coeff3*(-1 + mvc1*rhoc1*ur)));
dcoeff6 = h*dcoeff6;
