function [c, f, s, g0, pt] = flory_hugg_pde(x, t, u, dudx, a, b, e, u0,...
                                    e_g0, u_g0, v, j, mode, par_n)
% Solve with full ternary model. Analytical derivatives.
% pt ... phi_tot
% gra_a ... analytical gradient of phi_tot

pt = Ternary_model.phi_tot(x, a, b, e, u0, v*t);
gra_a = Ternary_model.gradient_analytical(x, a, b, e, v*t);
s = 0;
c = 1;
if strcmp(mode, 'Constituent')
    g0 = Ternary_model.gamma0(x, a+t*v, b, e_g0, u_g0, v*t, e, u0);
    f = g0.*(1-pt)./pt.*(pt.*dudx-u.*gra_a);
elseif strcmp(mode, 'Const_mob')
    % Using Stefano's mobility ansatz (1-phi_tot+a*phi_tot^2) to obtain a
    % certain ratio of D_in/D_out.
    f = (1-pt+e_g0*pt^2)./pt.*(pt.*dudx-u.*gra_a);
elseif strcmp(mode, 'Client')
    g0 = 0.5;
    chi_phi = -4.530864768482371;
    f = g0*(dudx+chi_phi*u*gra_a);
elseif strcmp(mode, 'Const_flux')
    g0 = Ternary_model.gamma0(x, a+t*v, b, e_g0, u_g0, v*t, e, u0);
    f = g0.*(1-pt)./pt.*(pt.*dudx-u.*gra_a)-j*u/pt;
elseif strcmp(mode, 'Mobility_ratio')
    f = u_g0.*(1-e_g0*pt).*(dudx-u.*gra_a./pt);
elseif strcmp(mode, 'Hammer_Ribo')
    % Use code from client (flux) and constituent (mobility)
    c = [1; 1];
    g01 = Ternary_model.gamma0(x, a+t*v, b, e_g0, u_g0, v*t, e, u0);
    g02 = Ternary_model.gamma0(x, par_n.a+t*par_n.v, par_n.b,...
                               par_n.e_g0, par_n.u_g0, ...
                               par_n.v*t, par_n.e, par_n.u0);

    f = [g01;g02].*(dudx+[par_n.chi_s; par_n.chi_p].*u*gra_a);
    reaction_rate = par_n.r * Ternary_model.phi_tot(x, a, b, ...
                                                       par_n.y, u0, v*t);
    F = reaction_rate*u(1);
    s = [-F; F];
end
end