function [f,p]=mcricc(beta,A1,A2,B,Q,R,W,equil)
m1=max(size(A1));
[row,m2]=size(A2);
Q11=Q(1:m1,1:m1);
Q12=Q(1:m1,m1+1:m1+m2);
Q21=Q(m1+1:m1+m2,1:m1);
Q22=Q(m1+1:m1+m2,m1+1:m1+m2);
W1=W(1:m1,:);
W2=W(m1+1:m1+m2,:);
p110=-.01*eye(m1);
p120=zeros(m1,m2);
p210=zeros(m2,m1);
p220=-.01*eye(m2);
dd=1;
it=1;
maxit=1000;
% check tolerance; for greater accuracy set it to 1e-10
while (dd>1e-10 & it<= maxit);
  f10=(2*R+beta*B'*(p110+p110')*B)\(2*W1'+beta*B'*(p110+p110')*A1);
  f20=(2*R+beta*B'*(p110+p110')*B)\(2*W2'+beta*B'*(p110+p110')*A2);
  f30=(2*R+beta*B'*(p110+p110')*B)\(beta*B'*(p120+p210'));
  eval(['g=',equil,'([f10 f20 f30]);'])
  G1=(eye(m2)+g*B*f30)\(g*(A1-B*f10));
  G2=(eye(m2)+g*B*f30)\(g*(A2-B*f20));
  AA1=A1-B*f10-B*f30*G1;
  AA2=A2-B*f20-B*f30*G2;
  F1=f10+f30*G1;
  F2=f20+f30*G2;
  p111=Q11+(F1'*R+2*W1)*F1+beta*(AA1'*p110*AA1+G1'*p210*AA1+AA1'*p120*G1+ ..
       G1'*p220*G1);
  p121=Q12+(F1'*R+2*W1)*F2+beta*(AA1'*p110*AA2+G1'*p210*AA2+AA1'*p120*G2+ ..
       G1'*p220*G2);
  p211=Q21+(F2'*R+2*W2)*F1+beta*(AA2'*p110*AA1+G2'*p210*AA1+AA2'*p120*G1+ ..
       G2'*p220*G1);
  p221=Q22+(F2'*R+2*W2)*F2+beta*(AA2'*p110*AA2+G2'*p210*AA2+AA2'*p120*G2+ ..
       G2'*p220*G2);
  f11=(2*R+beta*B'*(p111+p111')*B)\(2*W1'+beta*B'*(p111+p111')*A1);
  f21=(2*R+beta*B'*(p111+p111')*B)\(2*W2'+beta*B'*(p111+p111')*A2);
  f31=(2*R+beta*B'*(p111+p111')*B)\(beta*B'*(p121+p211'));
  dd=max(max(abs([f11 f21 f31]-[f10 f20 f30])./ ...
                max(ones([f10 f20 f30]),abs([f10 f20 f30])) ));
  it=it+1;
  p110=p111;
  p120=p121;
  p210=p211;
  p220=p221;
end;
f=[f11 f21 f31];
p=[p111 p121;
   p211 p221];
if it>=maxit; disp('WARNING: Iteration limit of 1000 reached in DRICC.M'); end;
