function [X, Y]=euler(x0, y0, xf, n)
    deff('yd=f(x,y)', 'yd=y-x^2+1');  
    h = (xf-x0)/n
    X = zeros(1,n+1)
    Y = X
    X(1) = x0
    Y(1) = y0
    xi = x0
    yi = y0
    for i=1:n   
        yi = yi + h*f(xi,yi)
        xi = xi+h
        Y(i+1) = yi
        X(i+1) = xi
    end
endfunction