function  jacob = vanvusse_Amatrix(x,par)
%definições das variáveis
Ca=  x(1);
Cb = x(2);
T  = x(3);
Tk = x(4);
%parâmetros externos
F  = par(1);
Qk = par(2);

%parâmetros comuns



%parâmetros cinéticos da eq Arrhenius
k0 = [1.287e12 1.287e12 9.043e9];
E_R = [9758.3 9758.3 8560.0];
k = k0.*exp(-E_R/T);
%delta H
DH = [-4.2 11 41.85];

V  = 10;
rho = 0.9342;
Cp = 3.01;
Kw = 4032;
Ar = 0.215;
T0 = 130 + 273.15;
Ca0 = 5.1;


df(1) = (F/V)*(Ca0 - Ca) - k(1)*Ca - k(3)*(Ca^2);
df(2) = -(F/V)*Cb +  k(1)*Ca - k(2)*Cb;
df(3) = (1/(rho*Cp))*(k(1)*Ca*DH(1) + k(2)*Cb*DH(2) + k(3)*Ca*Ca*DH(3)) +  (F/V)*(T0 - T) + ((Kw*Ar)/(rho*Cp*V))*(Tk - T);

%definições para o jacobiano
E1_R = E_R(1);
E2_R = E_R(2);
E3_R = E_R(3);
k10 = k0(1);
k20 = k0(2);
k30 = k0(3);
DH1 = DH(1);
DH2 = DH(2);
DH3 = DH(3);

mk = 5.0;   
cpk = 2.0; 

%jacobiano
%primeira linha
J(1,1) = - F/V - k10/exp(E1_R/T) - (2*Ca*k30)/exp(E3_R/T);
J(1,2) = 0;
J(1,3) = - (Ca^2*E3_R*k30)/(T^2*exp(E3_R/T)) - (Ca*E1_R*k10)/(T^2*exp(E1_R/T));
J(1,4) = 0;

%segunda linha
J(2,1) = k10/exp(E1_R/T);
J(2,2) = - F/V - k20/exp(E2_R/T);
J(2,3) = (Ca*E1_R*k10)/(T^2*exp(E1_R/T)) - (Cb*E2_R*k20)/(T^2*exp(E2_R/T));
J(2,4) = 0;

%terceira linha
J(3,1) = ((DH1*k10)/exp(E1_R/T) + (2*Ca*DH3*k30)/exp(E3_R/T))/(Cp*rho);
J(3,2) = (DH2*k20)/(Cp*rho*exp(E2_R/T));
J(3,3) = ((Ca*DH1*E1_R*k10)/(T^2*exp(E1_R/T)) + (Cb*DH2*E2_R*k20)/(T^2*exp(E2_R/T)) + (Ca^2*DH3*E3_R*k30)/(T^2*exp(E3_R/T)))/(Cp*rho) - F/V - (Ar*Kw)/(Cp*V*rho);
J(3,4) = (Ar*Kw)/(Cp*V*rho);

%quarta linha
J(4,1) = 0;
J(4,2) = 0;
J(4,3) = (Ar*Kw)/(cpk*mk);
J(4,4) = -(Ar*Kw)/(cpk*mk);


jacob = J;

 

