clear; close; clc;

% Define function
f = @(x) x.^3 - 5*x + 1;

% Initial bracket
a = 0; b = 2;
tol = 1e-4; maxit = 100;

for k = 1:maxit
    c = (a + b)/2;
    if abs(f(c)) < tol
        break
    elseif f(a)*f(c) < 0
        b = c;
    else
        a = c;
    end
end

root = c;
fprintf('Root = %.4f found in %d iterations\n', root, k)

clear; close; clc;

f = @(p) 20*p.^(-1.5) - (2 + 3*p);
a = 0.5; b = 5;
tol = 1e-4; maxit = 100;


for k = 1:maxit
    c = (a + b)/2;
    if abs(f(c)) < tol
        break
    elseif f(a)*f(c) < 0
        b = c;
    else
        a = c;
    end
end

root = c;
fprintf('Root = %.4f found in %d iterations\n', root, k)

clear; close; clc;

% f and derivative
f  = @(x) x.^3 - 5*x + 1;
fp = @(x) 3*x.^2 - 5;

x  = 1.5;                 % initial guess
tol = 1e-8; maxit = 50;

for it = 1:maxit
    fx  = f(x);
    fpx = fp(x);
    if abs(fx) < tol, break; end
    if abs(fpx) < 1e-12, warning('Derivative near zero'); break; end
    x = x - fx/fpx;       % update (alpha=1)
end

fprintf('Root ~ %.10f (it=%d, |f|=%.1e)\n', x, it, abs(f(x)));


alpha = 0.5;              % step size in (0,1] to stabilize
x = 1.5; tol = 1e-8; maxit = 50;

for it = 1:maxit
    fx  = f(x); fpx = fp(x);
    if abs(fx) < tol, break; end
    if abs(fpx) < 1e-12, warning('flat slope'); break; end
    xnew = x - alpha * fx / fpx;
    if ~isfinite(xnew), warning('bad update'); break; end
    x = xnew;
end

fprintf('Root ~ %.10f (it=%d, |f|=%.1e)\n', x, it, abs(f(x)));

%%
clear; close; clc;

xstar = fzero(@(x) x.^3 - 5*x + 1, 1.5);

F = @(z) [ z(1)^2 + z(2) - 2;
            exp(z(1)) + z(2) - 3 ];
z0 = [1; 1];
zstar = fsolve(F, z0);

%%

clear; close; clc;

w=1; a=0.5; sigma=2; chi=3.7; phi=2; tau=0; g=0;
Z  = @(h) ((1-tau)*w*h + a-g).^(-sigma).*(1-tau)*w - chi*(1-h).^phi;
dZ = @(h) -sigma*((1-tau)*w*h + a - g).^(-sigma-1).*(1-tau)*w + chi*phi*(1-h).^(phi-1);

h = 0.5; alpha = 1; tol = 1e-8;
for it=1:60
    Zh = Z(h); if abs(Zh) < tol, break; end
    dZh = dZ(h); if abs(dZh) < 1e-10, break; end
    h = h - alpha * Zh/dZh;           % damped Newton
    h = min(max(h, 0+1e-8), 1-1e-8);  % project to (0,1)
end
fprintf('h* = %.6f, Z=%.1e, it=%d\n', h, Z(h), it);

