clear       // clear all previous variables

// This only works for polynomials! We'll need other techniques for finding derivatives with Scilab, which we'll see later

function newtonPoly(x0,tolerance,maxstep)
    // The function takes three inputs: the starting guess of the root, 
    // the required tolerance, and the max number of steps at which to terminate in case
    // the method does not converge before then (always a good idea)

    
    // Draw a picture of the function near the starting guess:
    clf
    x = [ x0-10 : 0.1 : x0+10 ]
    plot(x,horner(f,x))
    g = gca();  
    g.x_location = "origin"; 
    g.y_location = "origin"; 
    plot(x0,horner(f,x0),'s')
 
 
    // Now for the method itself:
    step = 1        // Keep track of which step we're on
 
    converge = %F   // Assume the process won't work, change later if it does
    
    
    
     
    while (step <= maxstep)     // Keep repeating, unless maximum steps reached   
 
        df = derivat(f)         // derivative of polynomial f(x)
        
        num = horner(f,x0)      // Evaluate polynomial at x0
        den = horner(df,x0)     // Evaluate derivative at x0
        
        p = x0 - num/den;       // next root approximation
        fp = horner(f,p);       // value of f at next root approximation    
        
         
        if(fp == 0 | abs(fp) < tolerance) then
            disp('Step '+string(step))
            disp('     x_i: '+string(x0) )
            disp('     Next guess, x_i+1: '+string(p))
            disp('     f(x_i+1): '+string(fp))

            converge = %T
            break
        else
            // no root found yet, print step info
            disp('Step '+string(step))
            disp('     x_i: '+string(x0) )
            disp('     Next guess, x_i+1: '+string(p))
            disp('     f(x_i+1): '+string(fp))
            // add a marker to the axis
            
            plot(p,horner(f,p),'d')
        end
         
        // increase step counter  
        step = step + 1
         
        // redefine starting pair:
        x0 = p;
        
    end

    // print summary
    if (converge) then
        disp('Approximate root: '+string(p)+' found in '+string(step)+' steps')
    else
        disp('Root not found! Try different starting points or adjust tolerance')
    end
     
endfunction


f = poly( [ -2 0 1 ], 'x', 'coeff')     // define polynomial with coefficients (REVERSE standard form)

newtonPoly( -5, 0.0001, 100 )
disp(f)      // print polynomial
