function [L, G, i] = myeigorder(M)

%This function does the same thing Chris House's eig_order function does
%but since his function hasn't been posted, we wrote our own.

%Input: 
%       M - a square matrix
%Output: 
%       L - The eigenvalues of M along the diagonal 
%           increasing in absolute value
%       G - The eigenvectors of M with the same ordering as L.
%       i - The number of stable eigenvalues, also the index number
%           for the diagonal element that is the last stable eigenvalue.

N = size(M);

%Decompose the M matrix using Spectral Decomp.
[V,D] = eig(M);
%Calculate the absolute value of the lambda matrix and put them in a vector
Dabs = abs(D);
Dabsv = diag(Dabs);
%Sort the vector keeping the original index location
[L0, index] = sort(Dabsv,1);
%Put all the eigenvalues and their signs into vectors
Dsign = D < 0;
Dsign = -1*Dsign;
Dsignv = diag(Dsign);
Dsignv = (Dsignv == 0) + Dsignv;
%Use the original location and sign to recreate the L matrix with the new
%order. the command diag(v) puts the vector v on the main diagonal of a
%matrix or extracts a vector from the diagonal of a matrix.
L = diag(L0.*Dsignv);
%Use the original index location to reorder the eigenvector matrix V.
G = V(:,index');
%Number of stable eigenvectors
i = sum((abs(diag(L))<1 ==1),1);