function [c, f, s, g0, pt] = flory_hugg_pde(x, t, u, dudx, a, b, e, u0,...
                                    e_g0, u_g0, v, j, mode)
% 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);
end
end