with(plots):

# Systeemparameters
L := 200:  # Lengte van het systeem
K := 99:  # Aantal ruimtelijke punten
j := 2:  # Modusnummer
omega_d := 1:  # Kenmerkende frequentie
beta := 1:  # Niet-lineariteit parameter
delta := 0.05:  # Dempingscoëfficiënt

# Ruimtelijke raster
h := L / (K + 1):  # Stapgrootte
n := Vector(1..K+2, k -> -L/2 + (k-1)*h):  # Ruimtelijke punten
N := K + 2:  # Aantal punten inclusief randen

# Geschaalde parameters
omegaDScaled := h * omega_d:
deltaScaled := h * delta:

# Tijdparameters
dt := 1:  # Tijdstap
tmax := 3000:  # Maximale tijd
tspan := [seq(0, tmax, dt)]:  # Tijdvector

# Definieer de differentiaalvergelijkingsfunctie (odefun)
odefun := proc(t, Y, N, h, omegaDScaled, deltaScaled, beta)
    local U, Udot, Uddot, dUdt, dUdotdt, dYdt, k;
    
    U := Vector(1..N, k -> Y[k]):  # Verplaatsingen U
    Udot := Vector(1..N, k -> Y[k+N]):  # Snelheden Udot
    
    Uddot := Vector(1..N):  # Tweede afgeleiden (Laplaciaan)
    
    # Laplacian (discrete second derivative) voor interne punten
    for k from 2 to N-1 do
        Uddot[k] := (U[k+1] - 2*U[k] + U[k-1]);
    end do;
    
    # Randvoorwaarden
    Uddot[1] := 0;
    Uddot[N] := 0;
    
    # Stelsel van vergelijkingen
    dUdt := Udot;
    dUdotdt := Uddot - deltaScaled * Udot + omegaDScaled^2 * (U - beta * U^~3);
    
    # Pack afgeleiden in één vector dYdt
    dYdt := Vector(2*N):
    for k from 1 to N do
        dYdt[k] := dUdt[k]:
        dYdt[k+N] := dUdotdt[k]:
    end do:

    
    return dYdt;
end proc:

###########################################################################################
# Begincondities voor initieel amplitude a = 2
#a_init := 2:
#U0_init := Vector(1..N, k -> a_init * sin((j * Pi * h * n[k]) / L)):  # Beginverplaatsing
#U0_init[1] := 0:  # Randvoorwaarde aan de linkerzijde
#U0_init[N] := 0:  # Randvoorwaarde aan de rechterzijde
#Udot0_init := Vector(1..N, k -> 0):  # Begin snelheid = 0

# Plot de beginverplaatsing U0_init
#plot(U0_init, n, color=red, linestyle=solid, title="t=0, a=2", labels=["x_n", "U_n"]);
# Begincondities voor initieel amplitude a = 2
a_init := 2:
U0_init := Vector(1..N, k -> a_init * sin((j * Pi * h * n[k]) / L)):  # Beginverplaatsing blijft hetzelfde
U0_init[1] := 0:  # Randvoorwaarde aan de linkerzijde
U0_init[N] := 0:  # Randvoorwaarde aan de rechterzijde
Udot0_init := Vector(1..N, k -> 0):  # Begin snelheid = 0

# Plot de beginverplaatsing U0_init zonder de fysische betekenis te veranderen
plot(n, U0_init, color=red, linestyle=solid, title="t=0, a=2 (ongewijzigde fysische betekenis)", labels=["x_n", "U_n"],
     view=[-100..100, -3..3]);


#############################################################################################
# Waarden van amplitude 'a' om over te itereren
a_values := [2, 1.95, 1.9, 1.85, 1.82]:

# Loop over elke waarde van 'a'
for i to numelems(a_values) do
    a := a_values[i]:
    
    # Begincondities voor verplaatsing en snelheid
    U0 := Vector(1..N, k -> a * sin((j * Pi * h * n[k]) / L)):
    U0[1] := 0:  # Randvoorwaarde aan de linkerzijde
    U0[N] := 0:  # Randvoorwaarde aan de rechterzijde
    Udot0 := Vector(1..N, k -> 0):  # Beginsnelheid = 0

    # Pack begincondities
    Y0 := <U0 | Udot0>:

    # Oplossen van de ODE's
    sol := dsolve({seq(diff(Y[k](t), t$2) = odefun(t, [seq(Y[j](t), j = 1..2*N)], N, h, omegaDScaled, deltaScaled, beta)[k], k = 1..2*N),
                   seq(Y  = Y0[k], k = 1..2*N)}, numeric, output=listprocedure):

    # Extract de oplossing voor verplaatsing
    Usol := eval([seq(sol(t)[k], k=1..N)], t=tmax):
    
    # Plot de eindverplaatsing
    plot(Usol, n, color=blue, linestyle=solid, title=sprintf("t=%d, a=%.2f", tmax, a), labels=["x_n", "U_n"]);
end do;

Error, Matrix index out of range