% construct mass and stiffness matrices 
% compare eigenvalues to those of the Laplacian

 N = input('enter N: ');
 h = 1/(N+1);

% construct the mass matrix
 M = (2*h/3)*eye(N)+(h/6)*(diag(ones(N-1,1),1)+diag(ones(N-1,1),-1));

% construct the stiffness matrix
 K = (2/h)*eye(N)-(1/h)*(diag(ones(N-1,1),1)+diag(ones(N-1,1),-1));

% inspect the eigenvalues and eigenvectors of A = inv(M)*K

 A = M\K;                % A = inv(M)*K
 [V,D] = eig(A);         % compute eigenvalues and eigenvectors
 ew = diag(D);           % extract eigenvalues
 [ew,indx] = sort(ew);   % sort eigenvalues
 V = V(:,indx);          % reorder eigenvectors

 figure(1), clf, divax
 plot([1:N],ew,'b.','markersize',24)
 xlabel('$k$','fontsize',18,'interpreter','latex')
 ylabel('$\lambda_k$','fontsize',18,'interpreter','latex')
 xlim([1 N])
 set(gca,'xtick',[1 8 16 24 32]) 

 th = pi*[1:N]'/(N+1);

 for k=[1:N]
    v = sqrt(2/(N+1))*sin(k*th);
    v = v/max(abs(v)); 
    v = v*sign(v(1));
    figure(2), clf, divax
    plot([1:N],v,'b.-','markersize',24,'linewidth',1.5)
    xlabel('$j$','fontsize',18,'interpreter','latex')
    ylabel('$({\bf v_k})_j$','fontsize',18,'interpreter','latex')
    title(sprintf('eigenvector ${\\bf v}_k$ for $k=%d$',k),'interpreter','latex','fontsize',18)
    xlim([1 N])
    ylim(1.1*[-1 1])
    pause
 end
