function [V,D]= lazar_num(level,k,d,v,v1)
 r=d^((k-1)/k);  % renormalization factor for the energy 

%%  Calculate Previous Energy Matrix    %%
% ---------------------------------------------------------------------
% Here we create the matrix corresponding to the difference quotient on 
% the level prior 

% Get set up info 
%[rays1,vals1,weights1,matches1]=LapBW(level-1,creal,cimag,'n');
numpoints1=k*d^(level-1); 
rays1=v1(:,1);
% Calculate the differences btw ray angles- the denominators
%Drays1=[rays1,0]-[rays1(numpoints1)-((2^k-1)*d^(level-1)),rays1];
%Drays1=Drays1(1:numpoints1);
Drays1=zeros(numpoints1,1);
Drays1(1)=abs(rays1(numpoints1)-rays1(1));
for j=2:(numpoints1)
    Drays1(j)=abs(rays1(j-1)-rays1(j));
end
clear rays1;
% Calculate Old Energy Matrix -- the 1st difference quotient with theta
% factor
OE1=diag(-1./Drays1,0)+diag(1./Drays1(2:numpoints1),1);
OE1(numpoints1,1)=1/Drays1(1);
clear Drays1
OE1=6*OE1;  % To normalize our matrix 

% the 2nd difference quotient
OE2=diag(-1*ones(1,numpoints1),0)+diag(ones(1,numpoints1-1),-1);
OE2(1,numpoints1)=1;

OE=r^(level-1)*OE1*OE2; % take product and normalize 
clear OE1;
clear OE2;
%% Now Examine the Level of Interest  %%
% -----------------------------------------------------------
% repeat the process above but on the new level, also modify the previous
% energy matrix to opperate at this higher level

% Get info for this level 
rays=v(:,1);
numpoints=k*d^level;    

weights=ones(1,numpoints )/numpoints;
% Raise the energy matrix from the previous level to opperate at the new
% higher level
Eupped1=zeros(numpoints,numpoints/d);
Eupped1(~mod(rays,d),:)= OE; 
Eupped=zeros(numpoints,numpoints);
Eupped(:,~mod(rays,d))=Eupped1; 

% Make Matches Map (converts from vectors of free variables to
% eigenfunctions on the bascilica by making the necessary identifications).
M =zeros(numpoints,numpoints/k);
for j=1:numpoints
      M(j,v(j,2))=1;
end 

% Calculate the differences btw ray angles
%Drays=[rays,0]-[rays(numpoints)-((2^k-1)*d^level),rays];
%Drays=Drays(1:numpoints);
Drays=zeros(numpoints,1);
Drays(1)=abs(rays(numpoints)-rays(1));
for j=2:(numpoints)
    Drays(j)=abs(rays(j-1)-rays(j));
end
% Calculate Energy Matrix
E1=diag(-1./Drays,0)+diag(1./Drays(2:numpoints),1);
E1(numpoints,1)=1/Drays(1);
E1=6*E1;  % to get the normalization correct ?????
clear Drays;
clear rays;
E2=diag(-1*ones(1,numpoints),0)+diag(ones(1,numpoints-1),-1);
E2(1,numpoints)=1;

E=(1/2)*M'*(r^level*E1*E2+Eupped)*M;   %Use average of Eupped and E.  
clear E1;
clear E2;
clear Eupped;
% Calculate 'Integral' wrt measure matrix
G=M'*diag(weights)*M;
clear M;
% Laplacian Matrix
L=inv(G)*E;
clear G;
clear E;
% Calculate and order eigenvalues
[V,D]=eig(L);
clear L;
%e=vpa(e);
end
