Y = [0.12 -0.09 0.06 -0.05 0.032]';
T = [0.8 1.7 2.5 3.3 4.1]';

F = @(w,b) Y + 0.17 * exp(-b*T).*(cos(w*T) + b/w*sin(w*T));


update = 1;
tol = 1e-4;

w = 2*pi/(T(3)-T(1));
b = w / (2*pi) * (log(Y(1)) - log(Y(3)));
c = [w,b]';
s = 0.00001;

while norm(update, inf) > tol
    w = c(1); %Verbose
    b = c(2);
    ws = w + s;
    bs = b + s;
    Fw = F(ws,b);
    Fb = F(w,bs);
    Fval = F(w,b);
    DF = [(Fw-Fval)/s (Fb-Fval)/s];
    update = -DF\Fval;
    c = c + update

end


% w =3.7721
% b = 0.3903

w = c(1);
b = c(2);

tplot = linspace(0,6,1000);
Fplot = -0.17 * exp(-b*tplot).*(cos(w*tplot) + b/w*sin(w*tplot));
plot(T,Y,'o',tplot,Fplot,'m-')
