function [f_nm1,f_np1,phi_x] = VPFP_nonlinear(f_n,z)
global delv; global Nt; global X; global V;
%mesh(V, X, f_n(:,:,1))
[Nx,~,~] = size(f_n);[Nz, dim] = size(z);
h = fcn_h(z);
rho_n = reshape(sum(f_n,2)*delv,Nx,Nz);
for n = 1:Nt
    vdx_fn = vdx_fn_perd_1st(f_n);
    vdx_fn =  MultiD_times(vdx_fn,sig_1(z));
    rho_np1 = rho_1st(vdx_fn, rho_n);
    phi_x = Poisson(rho_np1 - h);
    [sqrt_M, sqrt_M_plus] = M_np1(phi_x);
    B = b_depz(vdx_fn,f_n,sqrt_M);
    B = MultiD_times(B,ones(Nz, dim)./sig_3(z));
    A_diag = A_VPFP_linear_diag(sqrt_M, sqrt_M_plus,z);
    f_np1 = f_1st_depz(B,A_diag,sqrt_M);
    f_n = f_np1;
    %rho_n = reshape(sum(f_np1,2)*delv,Nx,Nz);
    if n == Nt-1
        f_nm1 = f_np1;
    end
    %mesh(V, X, f_np1(:,:,1))
end
1;
