Here is my code in C for plotting the amplitude versus w/w0 for the Duffing equation :

#include #include #include

// Duffing Equation double f(double t, double x, double y, double gamma, double w0, double f0, double w, double beta) { return y; }

double g(double t, double x, double y, double gamma, double w0, double f0, double w, double beta) { return -2 * gamma * y - 2 * w0 * x - beta * x * x * x + f0 * cos(w * t); }

// RK4 Loop void rk4(double *t, double *x, double *y, double a, double b, double h, int N, double gamma, double w0, double f0, double w, double beta) { for (int k = 0; k < N - 1; k++) { double A1x = f(t[k], x[k], y[k], gamma, w0, f0, w, beta); double A1y = g(t[k], x[k], y[k], gamma, w0, f0, w, beta);

double A2x = f(t[k] + h / 2.0, x[k] + h / 2.0 * A1x, y[k] + h / 2.0 * A1y, gamma, w0, f0, w, beta); double A2y = g(t[k] + h / 2.0, x[k] + h / 2.0 * A1x, y[k] + h / 2.0 * A1y, gamma, w0, f0, w, beta);

double A3x = f(t[k] + h / 2.0, x[k] + h / 2.0 * A2x, y[k] + h / 2.0 * A2y, gamma, w0, f0, w, beta); double A3y = g(t[k] + h / 2.0, x[k] + h / 2.0 * A2x, y[k] + h / 2.0 * A2y, gamma, w0, f0, w, beta);

double A4x = f(t[k] + h, x[k] + h * A3x, y[k] + h * A3y, gamma, w0, f0, w, beta); double A4y = g(t[k] + h, x[k] + h * A3x, y[k] + h * A3y, gamma, w0, f0, w, beta);

x[k + 1] = x[k] + (h / 6.0) * (A1x + 2.0 * A2x + 2.0 * A3x + A4x); y[k + 1] = y[k] + (h / 6.0) * (A1y + 2.0 * A2y + 2.0 * A3y + A4y); t[k + 1] = t[k] + h; } }

double calculate_max(double *x, int N) { double max = -1e5; int w = 1000; for (int i = N - w; i < N; i++) x[i] > max ? max = x[i] : 1;

return max; }

// Main int main() { // Set integration interval double a = 0.0; double b = 1000.0;

// Set time step double h = 0.01;

// Differential equation parameters double gamma = 0.005; double w0 = 1.0; double f0 = 0.2; double beta = 1.5;

// Array dimensions int N = (int)((b - a) / h) + 1.0;

// Define arrays double *t = (double *)calloc(N, sizeof(double)); double *x = (double *)calloc(N, sizeof(double)); double *y = (double *)calloc(N, sizeof(double));

// Frequency range double w_ini = 0.1; double w_fin = 5.0; double w_pas = 0.001;

// Loop over different values of w for (double w = w_ini; w <= w_fin; w += w_pas) { // Initial conditions t[0] = a; x[0] = 1.0; y[0] = 0.0;

rk4(t, x, y, a, b, h, N, gamma, w0, f0, w, beta); double amplitude = calculate_max(x, N);

// Print frequency and corresponding amplitude printf("%e\t%e ", w / w0, amplitude); }

// Free allocated memory free(y); free(x); free(t);

return 0; }

How can I modify it to obtain this ?

Thank you