%DOUBLEJ2.M

function cov=doublej2(a1,b1,a2,b2);

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;

cov=gamma1;
