function giaiPTNewton()
    syms x;
    disp('Giai phuong trinh f(x)=0 bang phuong phap Newton...');
    f = input('Nhap ham f(x) = ');
    fprintf('\n');
    disp('Nhap khoang [a b]:');
    a = input('a = ');
    b = input('b = ');
    fprintf('\n');
    slv = 1;
    % Kiem tra dieu kien hoi tu

    if subs(f,x,a)*subs(f,x,b)>=0
        fprintf('Khoang cach ly nghiem khong hop le');
        slv = 0;
    elseif isreal(solve(diff(f,x)))
        x_ = solve(diff(f,x));
        for i=1:length(x_)
            if isreal(x_(i)) && (a<=x_(i)) && (x_(i)<=b)
                fprintf('Dao ham cua ham so co diem bang 0');
                slv = 0;
            end;
        end;
    end;
    if slv
        % Tim x0 theo Fourier
        f1 = diff(f,x);
        f2 = diff(f,x,2);
        if subs(f,x,a)*subs(f2,x,a)>0
            x0 = a;
        else x0 = b;
        end;

        % Xay dung day lap
        n = input('Nhap so lan lap: n = ');
        xCurrent = x0;
        xBefore = 0;
        for i=1:n
            xBefore = xCurrent;
            xCurrent = xCurrent - subs(f,x,xCurrent)/subs(f1,x,xCurrent);
        end;
        xDraw = linspace(a-2,b+2);
        yDraw = subs(f,x,xDraw);
        fprintf('Nghiem x%d = %f\n',n,xCurrent);
        yBefore = subs(f,x,xBefore);
        plot(xDraw,yDraw,'b-');
        hold on;
        % Tim tiep tuyen tai diem x(n-1)
        dif = subs(f1,x,xBefore);
        xa = a-2;
        syms y;
        ya = solve(y-yBefore-dif*(xa-xBefore),y);
        xb = b+2;
        yb = solve(y-yBefore-dif*(xb-xBefore),y);
        plot([xa,xb],[ya,yb],'-');
        grid on;
        plot(xBefore,yBefore,'r*');
        plot(xCurrent,0,'g*');
    end;



end