function [sqrt_M, sqrt_M_plus] = M_np1(phi_x)
global V;global delv;
[Nv,~] = size(V);
[Nx,Nz] = size(phi_x);
V_plus = [V(1)-delv;V;V(Nv)+delv];

phi_x_plus = reshape(repmat(phi_x,[Nv+2,1]),Nx,Nv+2,Nz);
phi_x = reshape(repmat(phi_x,[Nv,1]),Nx,Nv,Nz);


sqrt_M = 1/sqrt(2*pi)*exp(-(repmat(V',Nx,1,Nz) + phi_x).^2/2);
sqrt_M_plus = 1/sqrt(2*pi)*exp(-(repmat(V_plus',Nx,1,Nz) + phi_x_plus).^2/2);
sqrt_M = sqrt(sqrt_M);
sqrt_M_plus = sqrt(sqrt_M_plus);
