@Hullzie16  I tried now to replicate it with a loop but this way failed too,

Download v4.mw

@Hullzie16 As I can understand I have to create a pointplot spatial points over var:= seq(x[i](t), i = 1 .. km) the, but still I have not figured out how to replicate it

Download dsolve_fix_3.mw

@Hullzie16 These plots are created by Matlab for my research. You can read this post and maybe this will help you to understand what I want to create Numerical Simulation of the Discrete Klein-Gordon Equation: A Study on Damped, Driven Nonlinear Wave Systems with Spatially Extended Initial Conditions. 

@Hullzie16 I tried some ideas and Maple created some useful results, but still, I can not figure out how to create the above results. Thank you very much for your effort!

Download dsolve_fix_3.mw

Oh ok!! @acer I will have it on my mind. I will not do that again

@Hullzie16 When you try to plot it, do you receive results?

Using the following code 

odeplot(sol, [[t, x[100](t)]], t = 0 .. 3000, numpoints = 3500, title = "Trajectory of x[100](t)", axes = boxed, gridlines=true, color = ["Blue"], size = [600, 200]);

but I am looking for the plot a=2, at t=3000

@Hullzie16 Thank you, I appreciate that. The link that is attached is not working

If the symbolic approach is too cumbersome or doesn’t yield a solution that respects the boundary conditions well, a numerical solution is a more robust approach, especially when dealing with boundary conditions at infinity. Here's how you could proceed:

You can attempt to solve the ODE as a Bessel equation and apply boundary conditions to determine constants C1​ and C2.

Instead of dealing directly with x=−∞, choose a very large negative number as an approximation for infinity, say x=−100

Numerically, you can solve the ODE with approximated boundary conditions near −∞ using dsolve(..., numeric) in Maple.

1.pdf

Download solution.mw

@janhardo yes this code work perfectly! https://www.mathworks.com/matlabcentral/discussions/general/865216-numerical-simulation-of-the-discrete-klein-gordon-equation-a-study-on-damped-driven-nonlinear-wave

@janhardo I tried to do it on my own without getting the results I wanted 

@janhardo The following is the code I use on Matlab. But from now on I want to do everything in Maple so I'm trying to get the same result 

% Parameters
L = 200;  % Length of the system
K = 99;  % Number of spatial points
j = 2;  % Mode number
omega_d = 1;  % Characteristic frequency
beta = 1;  % Nonlinearity parameter
delta = 0.05;  % Damping coefficient

% Spatial grid
h = L / (K + 1);
n = linspace(-L/2, L/2, K+2);  % Spatial points
N = length(n);
omegaDScaled = h * omega_d;
deltaScaled = h * delta;

% Time parameters
dt = 1; % Time step
tmax = 3000; % Maximum time
tspan = 0:dt:tmax; % Time vector

% Values of amplitude 'a' to iterate over
a_values = [2, 1.95, 1.9, 1.85, 1.82];  % Modify this array as needed

% Differential equation solver function
function dYdt = odefun(~, Y, N, h, omegaDScaled, deltaScaled, beta)
    U = Y(1:N);
    Udot = Y(N+1:end);
    Uddot = zeros(size(U));

    % Laplacian (discrete second derivative)
    for k = 2:N-1
        Uddot(k) = (U(k+1) - 2 * U(k) + U(k-1)) ;
    end

    % System of equations
    dUdt = Udot;
    dUdotdt = Uddot - deltaScaled * Udot + omegaDScaled^2 * (U - beta * U.^3);

    % Pack derivatives
    dYdt = [dUdt; dUdotdt];
end

% Create a figure for subplots
figure;

% Initial plot
a_init = 2;  % Example initial amplitude for the initial condition plot
U0_init = a_init *  sin((j * pi * h * n) / L); % Initial displacement
U0_init(1) = 0; % Boundary condition at n = 0
U0_init(end) = 0; % Boundary condition at n = K+1
subplot(3, 2, 1);
plot(n, U0_init, 'r.-', 'LineWidth', 1.5, 'MarkerSize', 10); % Line and marker plot
xlabel('$x_n$', 'Interpreter', 'latex');
ylabel('$U_n$', 'Interpreter', 'latex');
title('$t=0$', 'Interpreter', 'latex');
set(gca, 'FontSize', 12, 'FontName', 'Times');
 xlim([-L/2 L/2]);
ylim([-3 3]);
grid on;

% Loop through each value of 'a' and generate the plot
for i = 1:length(a_values)
    a = a_values(i);

    % Initial conditions
    U0 = a * sin((j * pi * h * n) / L); % Initial displacement
    U0(1) = 0; % Boundary condition at n = 0
    U0(end) = 0; % Boundary condition at n = K+1
    Udot0 = zeros(size(U0)); % Initial velocity

    % Pack initial conditions
    Y0 = [U0, Udot0];

    % Solve ODE
    opts = odeset('RelTol', 1e-5, 'AbsTol', 1e-6);
    [t, Y] = ode45(@(t, Y) odefun(t, Y, N, h, omegaDScaled, deltaScaled, beta), tspan, Y0, opts);

    % Extract solutions
    U = Y(:, 1:N);
    Udot = Y(:, N+1:end);

    % Plot final displacement profile
    subplot(3, 2, i+1);
    plot(n, U(end,:), 'b.-', 'LineWidth', 1.5, 'MarkerSize', 10); % Line and marker plot
    xlabel('$x_n$', 'Interpreter', 'latex');
    ylabel('$U_n$', 'Interpreter', 'latex');
    title(['$t=3000$, $a=', num2str(a), '$'], 'Interpreter', 'latex');
    set(gca, 'FontSize', 12, 'FontName', 'Times');
      xlim([-L/2 L/2]);
ylim([-2 2]);
    grid on;
end

% Adjust layout
set(gcf, 'Position', [100, 100, 1200, 900]); % Adjust figure size as needed

I tried to solve the same problem another way; however, I cannot find a problem at the end of the code. Do you have any suggestions about this?

Download new_attempt.mw

@janhardo  The below worksheet has some changes that are needed but still, the plot is not the desired. Firstly, I want to be sure that the code gives me the desired results. In the initial post, there is what I want to create

Download v2.mw

@acer The animation is very useful indeed. Could I create a procedure that runs different values of alpha without changing alpha every time? 

For example alpha={2, 1.91, 1.9}

@C_R  I tried to plot the TSUCS2 Attractor with the command odeplot. However, I could not find a way to create the nice colorful output like with the DEplot above

Download TSUCS2.mw