Newton's Method on MATLAB for Stationary Solutions for the Non-linear Klein-Gordon Equation

Question

By interpreting the equation in this way, we can relate the dynamics described by the discrete Klein - Gordon equation to the behavior of DNA molecules within a biological system. I am trying to represent stationary points of the discrete Klein - Gordon equation. The results are as follows. Mathematica's are different from MATLAB's. What changes should I make to my Matlab code? Any suggestion?

% Parameters
numBases = 100; % Number of spatial points
omegaD = 0.2;  % Common parameter for the equation
% Preallocate the array for the function handles
equations = cell(numBases, 1);
% Initial guess for the solution
initialGuess = 0.01 * ones(numBases, 1);
% Parameter sets for kappa and beta
paramSets = [0.1, 0.05; 0.5, 0.05; 0.1, 0.2];
% Prepare figure for subplot
figure;
set(gcf, 'Position', [100, 100, 1200, 400]); % Set figure size

% Newton-Raphson method parameters
maxIterations = 1000;
tolerance = 1e-10;

for i = 1:size(paramSets, 1)
    kappa = paramSets(i, 1);
    beta = paramSets(i, 2);
    % Define the equations using a function
    for n = 2:numBases-1
        equations{n} = @(x) -kappa * (x(n+1) - 2 * x(n) + x(n-1)) - omegaD^2 * (x(n) - beta * x(n)^3);
    end
    % Boundary conditions with specified fixed values
    someFixedValue1 = 10; % Replace with actual value if needed
    someFixedValue2 = 10; % Replace with actual value if needed
    equations{1} = @(x) x(1) - someFixedValue1;
    equations{numBases} = @(x) x(numBases) - someFixedValue2;
    
    % Combine all equations into a single function
    F = @(x) cell2mat(cellfun(@(f) f(x), equations, 'UniformOutput', false));
    
    % Newton-Raphson iteration
    x_solution = initialGuess;
    for iter = 1:maxIterations
        % Evaluate function
        Fx = F(x_solution);
        
        % Evaluate Jacobian
        J = zeros(numBases, numBases);
        epsilon = 1e-6; % For numerical differentiation
        for j = 1:numBases
            x_temp = x_solution;
            x_temp(j) = x_temp(j) + epsilon;
            Fx_temp = F(x_temp);
            J(:, j) = (Fx_temp - Fx) / epsilon;
        end
        
        % Update solution
        delta = -J\Fx;
        x_solution = x_solution + delta;
        
        % Check convergence
        if norm(delta, inf) < tolerance
            fprintf('Convergence achieved after %d iterations.\n', iter);
            break;
        end
    end
    
    if iter == maxIterations
        error('Maximum iterations reached without convergence.');
    end
    
    % Plot the solution in a subplot
    subplot(1, 3, i);
    plot(x_solution, 'o-', 'LineWidth', 2);
    grid on;
    xlabel('n', 'FontSize', 12);
    ylabel('x[n]', 'FontSize', 12);
    title(sprintf('\\kappa = %.2f, \\beta = %.2f', kappa, beta), 'FontSize', 14);
end

% Improve overall aesthetics
sgtitle('Stationary States for Different \kappa and \beta Values', 'FontSize', 16); % Super title for the figure

Mathematica plots

numBases = 100;  (*Number of base pairs,representing a segment of DNA*)
kappa = 1/
   10;    (*Elasticity constant,chosen to reflect the \
stiffness;arbitrary unit*)
omegaD = 
  2/10;   (*Frequency term,speculative and low to reflect slow \
dynamics*)
beta = 5/
   100;    (*Nonlinearity coefficient small to introduce mild \
nonlinearity*)

(*Initialize a list for the x_n variables*)
variables = Table[x[n], {n, 1, numBases}];

(*Define the equations based on the discretized version of the \
Klein-Gordon equation*)
equations = 
  Table[-kappa (x[n + 1] - 2 x[n] + x[n - 1]) == 
    omegaD^2 (x[n] - beta x[n]^3), {n, 2, numBases - 1}];

(*Add boundary conditions,for example,x[1] and x[numBases] fixed to 0 \
or any other value*)
boundaryConditions = {x[1] == 10, x[numBases] == 10};

(*Combine the equations with the boundary conditions*)
fullSystem = Join[equations, boundaryConditions];

(*Solve the system with an initial guess for each variable*)
initialGuess = 
  Table[{x[n], 1/100}, {n, 1, 
    numBases}]; (*or any small value close to expected solutions*)
solution = 
  FindRoot[fullSystem, initialGuess, 
   Method -> {"Newton", "StepControl" -> "TrustRegion"}];
(*Extract the solution for plotting or analysis*)
solvedValues = Table[x[n] /. solution, {n, 1, numBases}];

(*Optionally,you can plot the solution to visualize the stationary \
states*)
ListPlot[solvedValues, PlotStyle -> PointSize[Medium], Joined -> True,
  GridLines -> Automatic, PlotRange -> All, 
 AxesLabel -> {"n", "x[n]"}, LabelStyle -> {"Medium", Bold}]