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);