function [M,C] = basicsol(A,n1,n2,SigmaXX,x10,cutoff );
%%
%%-----------------------------------------------------------------------
%%
%%  Basicsol -- based on Soderlinds SimpRule but modified for solving an
%%  RBC model (Paul Söderlind, Paul.Soderlind@unisg.ch, Aug 2000, July 2007)
%%
%%  Purpose:   Solving RE model 
%%  The shocks have the covariance matrix Cov[e(t)] = SigmaXX.
%%
%%  Usage:    [M,C] = basicsol(A,n1,n2,SigmaXX,x10,cutoff);
%%
%%  n1 = number or predetermined variables
%%  n2 = number of nonpredetermined (forward looking) variables
%% 
%%  basic set up is E_t(Z_t+1) = A * Z_t + V_t+1
%%
%%  Input:     M,bet,n1,n2
%%             
%%             SigmaXX    n1xn1 covariance matrix of shocks
%%             x10        n1x1 vector, initial values for predetermined
%%                        variables
%%             cutoff     scalar, max modulus of backward looking variables
%%
%%
%%  Output:    M      -   n1xn1 matrix, x1(t+1) = M*x1(t) + e(t+1)
%%             C          n2xn1 matrix, x2(t) = C*x1(t)
%%
%%
%%-----------------------------------------------------------------------*/


  n = n1 + n2;

 [S,T,Qa,Z] = qz(eye(size(A)),A); %MatLab: I=Q'SZ' and A=Q'TZ'; Paul S:  I=QSZ' and A=QTZ',%but Q isn't used
 %[S,T,QQ,Z] = reorder(S,T,QQ,Z);  % reordering of generalized eigenvalues in ascending order
 logconA = abs(diag(T)) <= (abs(diag(S))*cutoff);   %selecting stable eigenvalues
 [S,T,Qa,Z] = ordqz(S,T,Qa,Z,logconA);
 logcon = abs(diag(T)) <= (abs(diag(S))*cutoff);   %1 for stable eigenvalue

if sum(logcon) < n1;
  warning('Too few stable roots: no stable solution');
  M = NaN; C = NaN;
  return;
elseif sum(logcon) > n1;
  warning('Too many stable roots: infinite number of stable solutions');
  M = NaN; C = NaN;
  return;
end;

 Stt = S(1:n1,1:n1);
 Zkt = Z(1:n1,1:n1);
 Zlt = Z(n1+1:n,1:n1);
 Ttt = T(1:n1,1:n1);

  if (cond(Zkt) > 1e+14) | (cond(Stt) > 1e+14);
    error('Zkt or Stt is singular');
  end;
  Zkt_1 = inv(Zkt);         %%inverting
  Stt_1 = inv(Stt);

  M = real(Zkt*Stt_1*Ttt*Zkt_1);       %%x1(t+1) = M*x1(t) + e(t+1)*/
  C = real(Zlt*Zkt_1);                 %%x2(t) = C*x1(t)*/

%%----------------------------------------------------------------------*/

