function [f_nm1, f_np1] = VPFP_linear(f_n, z)
global Nt; global V; global X;
phi_x = fcn_phi_x(z);
[sqrt_M, sqrt_M_plus] = M_np1(phi_x);
[Nz,dim] = size(z);
A_diag = VPFP_A_diag(sqrt_M, sqrt_M_plus,z);
for n = 1:Nt
    vdx_fn = vdx_fn_perd_1st(f_n);
    vdx_fn =  MultiD_times(vdx_fn,sig_1(z));
    B = VPFP_B(vdx_fn,f_n,sqrt_M);
    B = MultiD_times(B,ones(Nz, dim)./sig_3(z));
    f_np1 = f_1st(B,A_diag,sqrt_M);
    if n == Nt
        f_nm1 = f_n;
    end
    f_n = f_np1;
    %mesh(V,X,f_n(:,:,1));
end
1;