% EXS 3.8
clear all
format compact
k1 = 0.12; k2 = 0.22;
atot = 0.60; btot = 1.0;
% a = atot/(1+k1b+k2b^2) = (btot-b)/(k1 b + 2k2 b^2)
% atot(k1 b + 2k2 b^2)-(btot-b)(1+k1b+k2b^2) = 0
% - btot + b(atot k1 +1 -btot k1) +
%  b^2 (2 k2 atot -btot k2 + k1) +
%  b^3 k2
c = [k2,(2*k2*atot-btot*k2+k1),(atot*k1+1-btot*k1),-btot ];
z = roots(c);
i = find(imag(z)==0); % there is only one real root
x = real(z(i))
x0 = x;
%
% Analytical differentiation
%===================================
% dp/dz* dz/dk + dp/dk = 0
% dz/dk = -dp/dk / dp/dz;
dc      = polyder(c);
dpdz    = polyval(dc,x);
%
dcdk1   = [0, 1          , atot-btot, 0 ];
dcdk2   = [1, 2*atot-btot, 0        , 0 ];
dcdatot = [0, 2*k2       , k1       , 0 ];
dcdbtot = [0, -k2        ,-k1       ,-1 ];
% dp/dk:
dpdk1   = polyval(dcdk1  ,x);
dpdk2   = polyval(dcdk2  ,x);
dpdatot = polyval(dcdatot,x);
dpdbtot = polyval(dcdbtot,x);
% dz/dk:
dzdk1   = -dpdk1  /dpdz
dzdk2   = -dpdk2  /dpdz
dzdatot = -dpdatot/dpdz
dzdbtot = -dpdbtot/dpdz
%
% Experimental evaluation; change parameters by h
%===============================
h = 1e-4; % how large?
par0 = [k1 k2 atot btot];
npar = length(par0);
dxdp = zeros(4,1);
for k = 1:npar
    par = par0;
    par(k) = par(k)+h;
    k1   = par(1);
    k2   = par(2);
    atot = par(3);
    btot = par(4);
    c = [k2,(2*k2*atot-btot*k2+k1),(atot*k1+1-btot*k1),-btot ];
    z = roots(c);
    i = find(imag(z)==0);
    x = real(z(i));
    dxdp(k) = (x-x0)/h;
end
dxdp
