function [root, iter] = newton(func, x0, tol, max_iter)
    iter = 0;
    h = 1e-6; % Small step for numerical derivative
    while iter < max_iter
        iter = iter + 1;
        f_val = func(x0);
        f_prime = (func(x0 + h) - f_val) / h; % Numerical derivative
        x_new = x0 - f_val / f_prime;
        if abs(x_new - x0) < tol
            break;
        end
        x0 = x_new;
    end
    root = x0;
end