syms c1 c2 V alpha positive real
syms L1 L2 wT wD kk real 

q = sym('q(c1,c2)'); R1 = 1 ; R2 = sym('R2(c1,c2)');

%kk = 1/(c0 * V0)
V = 1 / (c1*kk); %dVdt = (alpha-1)*q;
dc1dt = q*c1/V*(R1-alpha);
dc2dt = q*c2/V*(R2-alpha);

Ha = wD*q - kk*(L1*q*c1*c1+L2*q*c1*c2);
H0 = wT + kk*L1*q*c1*c1 + kk*L2*q*c1*c2*R2;
H = Ha * alpha + H0;
 
dL1dt = diff(-H, c1); dL2dt = diff(-H, c2);

dqdt = diff(q,c1)*dc1dt+diff(q,c2)*dc2dt;
dR2dt = diff(R2,c1)*dc1dt+diff(R2,c2)*dc2dt;

dH0dt = dL1dt*c1^2*kk*q + L1*2*c1*dc1dt*kk*q + L1*c1^2*kk*dqdt;
dH0dt = dH0dt + dL2dt*q*c1*c2*kk*R2 + L2*dqdt*c1*c2*kk*R2;
dH0dt = dH0dt + L2*q*dc1dt*c2*kk*R2 + L2*q*c1*dc2dt*kk*R2 + L2*q*c1*c2*kk*dR2dt;
dHadt = wD*dqdt - kk*(dL1dt*c1^2*q +L1*2*c1*dc1dt*q+L1*c1^2*dqdt);
dHadt = dHadt - kk*(dL2dt*c1*c2*q +L2*dc1dt*c2*q+L2*c1*dc2dt*q+L2*c1*c2*dqdt);

dH0dt_coef = coeffs(dH0dt,alpha);
dH0dt0 = dH0dt_coef(1); dH0dta = dH0dt_coef(2); 

dHadt_coef = coeffs(dHadt,alpha);
dHadt0 = dHadt_coef(1); dHadta = dHadt_coef(2);

%H0 = wT + kk*L1*q*c1*c1 + kk*L2*q*c1*c2*R2;
tmp = coeffs(H0,L2); b13 = tmp(2)/c2; % kk*q*c1*R2;
tmp = coeffs(tmp(1),L1); b11 = tmp(1); b12 = tmp(2)/c1;

%Ha = wD*q - kk*(L1*q*c1*c1+L2*q*c1*c2);
tmp = coeffs(Ha,L2); b23 = tmp(2)/c2; 
tmp = coeffs(tmp(1),L1); b21 = tmp(1); b22 = tmp(2)/c1;

%dH0dta
tmp = coeffs(factor(dH0dta/(c1*kk)),L2);
b33 = simplify(tmp(2)/c2);  tmp = coeffs(tmp(1),L1);
b31 = tmp(1); b32 = tmp(2)/c1;

b = [b11 b12 b13;b21 b22 b23;b31 b32 b33]; db = factor(det(b));

tmp = coeffs(db, wT); sT = tmp(2)
tmp = coeffs(db, wD); sD = tmp(2)
