
function V=doublej2(a1,b1,a2,b2);
%  DOUBLEJ2.M
%  function V=doublej2(a1,b1,a2,b2);
%
%  Computes the matrix V defined by the infinite sum
%
%      V = SUM   [a1^j]*[b1*b2]*[a2^j]
%
%  where the sum goes from j = 0 to j = infinity.
%  The following doubling algorithm is used.
%  Starting from V(0) = b1*b2, iterate to convergence of V(j) on
%
%       a1(j) = a1(j-1)*a1(j-1)
%       a2(j) = a2(j-1)*a2(j-1)
%
%       V(j) = V(j-1) + a1(j-1)*V(j-1)*a2(j-1)
%
%  The limit point is returned in V.
alpha01=a1;
alpha02=a2;

gamma0=b1*b2';

diff=5;

while diff > 1e-15;

  alpha11=alpha01*alpha01;
  alpha12=alpha02*alpha02;

  gamma1=gamma0+alpha01*gamma0*alpha02';
  diff=max(max(abs(gamma1-gamma0)));
  gamma0=gamma1;
  alpha01=alpha11;
  alpha02=alpha12;

end;

V=gamma1;
