[Ps,Pu,L,V]=proj(A)

function [Ps,Pu,L,V]=proj(A)
% calculates the projection operators, the eigenvalues 
% and the eigenvectors of the Jacobian of A

[v,d]=eig(A);
[L,ind]=sort(diag(d));% the eigenvalues
V=v(:,ind);% the eigenvectors
           % V(:,1),V(:,2) basis in E_s, V(:,3) basis in E_u

Ps=V*diag([1 1 0])/V;Pu=V*diag([0 0 1])/V;