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

  a = s;

  n = length(d);

  for k = 1:maxIter

  maxiter = 0;

  ++maxiter
  daux = d .- h;
  saux = zeros(1, length(s));

   saux(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);
    saux(i) = saux(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) - saux(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)*saux(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