%% Solving example 1 (equilibrium reaction)
%
% First, rearrange the equation in $x$ into a polynomial of the form
% $f(x) = c_1 x^n + c_2 x^{n-1} + \cdots + c_{n+1} = 0$
%
% Namely,
% $-4Kx^3+4K(B_{ini}+A_{ini})x^2-(4KA_{ini}B_{ini}+KA_{ini}^2+1)x+KA_{ini}^2B_{ini}-C_{ini} = 0$
%
% Then set parameter values and define external or anonymuos function.

K=0.016;
Aini=42;
Bini=28;
Cini=4;

equreac = @(x) (-4*K*x.^3+4*K*(Bini+Aini)*x.^2-(4*K*Aini*Bini+K*Aini^2+1)*x+K*Aini^2*Bini-Cini);

%%
% Plot function.

x = linspace(0,42);
plot(x,equreac(x))
grid

%%
% Calculate coefficients of the polynomial and use roots to solve this.

C = [-4*K 4*K*(Bini+Aini) -4*K*Aini*Bini-K*Aini^2-1 K*Aini^2*Bini-Cini];

x = roots(C)

%%
% We are only interested in the real solution (the imaginary solutions are
% unphysical).
%
% Test the solution by substituting into the equilibrium equation.

x = x(3);

(Cini+x)/(Aini-2*x)^2/(Bini-x)

%%
% We could also have used fzero.

x = fzero(equreac,0)




