% pcaixavar.m
clear all
nm=4; % A dimensão da base será nm+1
ni=0;np=1; %todos os estados
%ni=0;np=2; %somente par
%ni=1;np=2; %somente impar
i=0;j=0;
% Construção das matrizes H e S
for n=ni:np:nm,i=i+1;for m=ni:np:nm;j=j+1;
    if mod((m+n),2)==0 
            h(i,j)= -8.0*(1.0-m-n-2*m*n)/((m+n+3.0)*(m+n+1.0)*(m+n-1.0));
            s(i,j) = 2.0/(n+m+5.0) - 4.0/(n+m+3.0) + 2.0/(n+m+1.0);
    end,end,j=0;end
[vet,val]=eig(h,s);vv=diag(val);[ma,na]=size(vv);
n=ni+1:np:nm+1;en=n.^2*pi^2/4;res=[n' vv en'];

figure(1)
x=-1:.1:1;
phic=0.99999*(1-x.^2)-0.23343*x.^2.*(1-x.^2)+0.018959*x.^4.*(1-x.^2); %com mn=4
n1=trapz(x,phic.^2),n1=sqrt(n1),phic=phic/n1
plot(x,phic,'o'),hold on,plot(x,cos(pi*x/2))

figure(2)
x=-1:.1:1;
phic2=3.07*x.*(1-x.^2)-1.5876*x.^3.*(1-x.^2)
n2=trapz(x,phic2.^2),n2=sqrt(n2),phic2=phic2/n2
plot(x,phic2,'o'),hold on,plot(x,sin(pi*x))

