k1 = 0.12;
k2 = 0.22;
atot = 0.6;
btot = 1.0;

a = atot/2;
b = btot/2;
c = [a b]';
update = 1;

tol = 1e-4;


F = @(a,b) [a+k1*a*b+k2*a*b^2-atot;...
    b+k1*a*b+2*k2*a*b^2-btot];

DF = @(a,b) [1+k1*b+k2*b^2,k1*a+2*k2*a*b;k1*b+2*k2*b^2,1+k1*a+4*k2*a*b];

while norm(update,inf) > tol
    update = -DF(c(1),c(2))\F(c(1),c(2));
    c = c + update
    
end