diff --git a/matlab/AIM/SPEigQZ.m b/matlab/AIM/SPEigQZ.m
deleted file mode 100644
index eb5153dbcad9b1ed9795caf1bd2933745eec130d..0000000000000000000000000000000000000000
--- a/matlab/AIM/SPEigQZ.m
+++ /dev/null
@@ -1,30 +0,0 @@
-function [w,rts,lgroots] = SPEigQZ(a,uprbnd,rowsLeft) 
-%  [w,rts,lgroots] = SPEigQZ(a,uprbnd)
-%
-%  Compute the roots and the left eigenvectors of the companion
-%  matrix, sort the roots from large-to-small, and sort the
-%  eigenvectors conformably.  Map the eigenvectors into the real
-%  domain. Count the roots bigger than uprbnd.
-
-
-eigOpt.disp=0;
-[w,d]   = eig(full(a'));
-rts     = diag(d);
-mag     = abs(rts);
-[mag,k] = sort(-mag);
-rts     = rts(k);
-
-ws=SPSparse(w);
-ws       = ws(:,k);
-
-%  Given a complex conjugate pair of vectors W = [w1,w2], there is a
-%  nonsingular matrix D such that W*D = real(W) + imag(W).  That is to
-%  say, W and real(W)+imag(W) span the same subspace, which is all
-%  that aim cares about. 
-
-ws = real(ws) + imag(ws);
-
-lgroots = sum(abs(rts) > uprbnd);
-
-w=full(ws);
-
diff --git a/matlab/AIM/SPQZ.m b/matlab/AIM/SPQZ.m
deleted file mode 100644
index 2a0ee66b887358f5f361ce9511542d9d614a76cd..0000000000000000000000000000000000000000
--- a/matlab/AIM/SPQZ.m
+++ /dev/null
@@ -1,91 +0,0 @@
-function  [q,rts,ia,nexact,nnumeric,lgroots,aimcode] = ...
-                       SPQZ(aa,bb)
-%  [q,rts,ia,nexact,nnumeric,lgroots,aimcode] = ...
-%                       SPQZ(aa,bb)
-%
-%  Solve a linear perfect foresight model using the matlab eig
-%  function to find the invariant subspace associated with the big
-%  roots.  This procedure will fail if the companion matrix is
-%  defective and does not have a linearly independent set of
-%  eigenvectors associated with the big roots.
-% 
-%  Input arguments:
-% 
-%    h         Structural coefficient matrix (neq,neq*(nlag+1+nlead)).
-%    neq       Number of equations.
-%    nlag      Number of lags.
-%    nlead     Number of leads.
-%    condn     Zero tolerance used as a condition number test
-%              by numeric_shift and reduced_form.
-%    uprbnd    Inclusive upper bound for the modulus of roots
-%              allowed in the reduced form.
-% 
-%  Output arguments:
-% 
-%    b         Reduced form coefficient matrix (neq,neq*nlag).
-%    rts       Roots returned by eig.
-%    ia        Dimension of companion matrix (number of non-trivial
-%              elements in rts).
-%    nexact    Number of exact shiftrights.
-%    nnumeric  Number of numeric shiftrights.
-%    lgroots   Number of roots greater in modulus than uprbnd.
-%    aimcode     Return code: see function aimerr.
-[neq,aCols]=size(aa);
-[bRows,bCols]=size(bb);
-if(~and(neq ==aCols,and(neq== bRows,neq ==bCols)))
-error(SPQZ:'aa, bb must be square and same dim')
-end
-nlag=0
-nlead=1
-b=[];
-rts=[];
-ia=[];
-nexact=[];
-nnumeric=[];
-lgroots=[];
-aimcode=[];
-
-condn=1.0e-8
-uprbnd=1.0e16
-
-% Initialization.
-nexact   = 0;
-nnumeric = 0;
-lgroots  = 0;
-iq       = 0;
-aimcode    = 0;
-b=0;
-qrows = neq*nlead;
-qcols = neq*(nlag+nlead);
-bcols = neq*nlag;
-q        = zeros(qrows,qcols);
-rts      = zeros(qcols,1);
-
-% Compute the auxiliary initial conditions and store them in q.
-h=[aa bb];
-
-
-[h,q,iq,nexact] = SPExact_shift(h,q,iq,qrows,qcols,neq);
-   if (iq>qrows) 
-      aimcode = 61;
-      return;
-   end
-
-[h,q,iq,nnumeric] = SPNumeric_shift(h,q,iq,qrows,qcols,neq,condn);
-   if (iq>qrows) 
-      aimcode = 62;
-      return;
-   end
-
-%  Build the companion matrix.  Compute the stability conditions, and
-%  combine them with the auxiliary initial conditions in q.  
-
-[a,ia,js] = SPBuild_a(h,qcols,neq);
-
-if (ia ~= 0)
-   [w,rts,lgroots] = SPEigQZ(a,uprbnd,min(length(ia),qrows-iq+1));
-%enough in w to fill q
-   q = SPCopy_w(q,w,js,iq,qrows);
-end
-
-
diff --git a/matlab/AIM/SPSparse.m b/matlab/AIM/SPSparse.m
deleted file mode 100644
index 294b19e25cf90626c4a7074c18518bd8fc60c771..0000000000000000000000000000000000000000
--- a/matlab/AIM/SPSparse.m
+++ /dev/null
@@ -1,2 +0,0 @@
-function res = SPSparse(aMat)
-res=sparse(aMat);