function [valp vecp] = PutereInv(d, s, h, y, maxIter, tol)

  a = s

  n = length(d)

  for k = 1:maxIter

  daux = d .- h

   s(1) = s(1)/daux(1)
   y(1) = y(1)/daux(1)


   for i = 2:n-1
    t = daux(i) - a(i-1)*s(i-1)
    s(i) = s(i)/t
    y(i) = (y(i) - y(i-1)*a(i-1))/t
   endfor

   y(n) = (y(n) - y(n-1)*a(n-1))/t;

   z(n) = y(n);

   for i = n-1:-1:1
    z(i) = y(i) - s(i)*z(i+1);
   endfor

   vecp = z/norm(z)

   valp = 0;
   for i = 1:n
    valp += (vecp(i)^2)*d(i)
   endfor

   for i = 1:n-1
    valp += 2*(vecp(i)*vecp(i+1)*s(i))
   endfor

   temp = zeros(n, 1);

   temp(1) = d(1)*vecp(1) + s(1)*vecp(2)
   temp(n) = s(n-1)*vecp(n-1) + d(n)*vecp(n)

   for i = 2:n-1
    temp(i) = s(i-1)*vecp(i-1) + vecp(i)*d(i) + s(i)*vecp(i+1)
   endfor

   if(norm(temp - valp*vecp) < tol)
    break;
   endif

   h = valp;

  endfor

endfunction