clear, clc;

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

F=@(atot,btot,k1,k2,a,b) [a+k1*a*b+k2*a*b^2-atot; b+k1*a*b+2*k2*a*b^2-btot];
J=@(k1,k2,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];

X0=[atot;btot];
while norm(F(atot,btot,k1,k2,X0(1),X0(2)))>1e-4
    f=F(atot,btot,k1,k2,X0(1),X0(2));
    j=J(k1,k2, X0(1),X0(2));
    delta=-j\f;
    X0=X0+delta;
end
X0
F(atot,btot,k1,k2,X0(1),X0(2))
