Transition point discontinuity (cliff/spike) in plots of solved system of non linear equations.

I have the following code to plot three graphs after solving system of non linear equations (Using the help of chat gpt). The axial coordinate is split in to two regions, 1. clean region 0≤x≤x_cr. and cake region x>x_cr. The system of equations is solved after certain x value x_cr while for the x ≤ x_cr, there is initial boundary condition imposed as follows J(x)=0, e(x)=0, and ε(x)=ε0.

The problem

When plotting J(x), e(x), and ε(x), I almost always observe a sharp edge / spike / cliff immediately after x_cr. This becomes more pronounced when changing parameters such as TMP and V. Physically, these profiles are expected to be smooth, so I believe the artifact is numerical, not physical.

What I am asking for help

I would greatly appreciate advice on how to reformulate the numerical treatment of the transition at x_cr to remove this discontinuity,

in particular:

• How should G be initialized and accumulated across a piecewise-defined domain?

• Should the first cake point be solved differently (e.g., special initialization at x_cr)?

• Is my use of trapz inside the marching loop appropriate, or should G_xi be advanced incrementally?

• Is there a better way to enforce continuity at the clean–cake interface?

Specifically, I am looking for guidance on correcting how G_xi and J(x) are computed and initialized near x_cr to avoid the observed spike.

The figure below shows one of the three graph with spike (which I need to remove by adjusting the code).

I pasted part of the code which contains the solve block and G_xi for trapizodal method for your reference.

I really appreciate your help already.

Kind regards,

Tadele

%% ---------------- Define regions clean vs cake ----------------
idx_cr = find(x > x_cr, 1, 'first');   % STRICTLY greater than x_cr
if isempty(idx_cr)
    error('x_cr is beyond L. Increase L or check parameters.');
end

% Clean region: up to idx_cr (inclusive)
J(1:idx_cr)   = J0;
e(1:idx_cr)   = 0;
eps(1:idx_cr) = eps_clean;

% Start solving only AFTER idx_cr
start_i = idx_cr + 1;
if start_i > Nx
    warning('No cake region points after x_cr for this discretization.');
end

%% ---------------- fsolve options ----------------
e_guess_prev = 1e-8;

opts = optimoptions('fsolve', ...
    'Display','off', ...
    'TolFun',1e-12, ...
    'TolX',1e-12, ...
    'MaxFunctionEvaluations', 2000, ...
    'MaxIterations', 2000);

%% ---------------- Cake region solve: i = start_i .. Nx ----------------
for i = start_i:Nx
    xi = x(i);

    % G(xi) = ∫0^xi J(s) ds using upstream known J values
    G_xi = trapz(x(1:i-1), J(1:i-1));
    if G_xi <= 0
        G_xi = J0 * xi;
    end

    e_guess = max(min(e_guess_prev, 0.9*R), 1e-9);

    % Residual uses bidisperse a_eff4 in Eq.5
    F = @(e_local) flux_diff_e_romero_davis_bidisperse( ...
        e_local, G_xi, ...
        Qcr0, phi0, mu0, a_eff4, ...
        TMP, Rm, rho, shape, d_c, R, ...
        tau_fun, eps_pref, eps_exp);

    [e_sol, ~, exitflag] = fsolve(F, e_guess, opts);
 

    e(i) = e_sol;
    e_guess_prev = e_sol;

    % Eq. 5 (your form), but with a_eff4
    tau_x = tau_fun(e_sol);
    J_eq5 = sqrt( Qcr0 * tau_x^3 * a_eff4 / (phi0 * mu0^3 * G_xi) );
    J(i)  = J_eq5;
    %++++++++++++++++++++++++++++++++++
    % Eq. 9: epsilon (same as your code)
     eps(i) = eps_pref;
    
end

Answer 1

[Matt J]

I find it strange that you are using fsolve for a 1-variable problem, rather than fzero. You don't show the objective function, so we don't even know if F() satisfies the continuity and differentiability assumptions of fsolve.

As for the discontinuity, you could try adding bounds on the search, to keep the solution close to the previous one. That would also be very easy with fzero, e.g.,

e_sol = fsolve(F, e_guess, e_guess+1e-3*[-1,+1]);

Answer 2

[Image Analyst]

You could post-process your curve data to remove the spike. Perhaps findchangepts will be useful to you.