function A_diag = VPFP_A_diag(M, M_plus,z)
global eps; global delv; global delt;
[Nx,Nv,Nz] = size(M);
%M = ones(Nx,Nv,Nz)./M;
%M_plus = ones(Nx,Nv+2,Nz)./M_plus; 
sig = reshape(repmat(sig_3(z),[1,Nx*Nv])',Nx,Nv,Nz);
A_diag = -(M_plus(:,3:Nv+2,:)./M + M_plus(:,1:Nv,:)./M + eps*delv^2/delt*ones(Nx,Nv,Nz)./sig);
