Here is my code to plot in c the amplitude versus w/w0 for the Duffing equation without the x^3 dependecy :

#include #include #include

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

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

void rk4(double *t, double *x, double *y, double a, double b, double h, int N, double gamma, double w0, double f0, double w) { for (int k = 0; k

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

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

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

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; } }

int main() { // Integration interval double a = 0.0; double b = 500.0;

// Time steps double h = 1.0;

// Differential equation parametres double gamma = 0.1; double w0 = 1.0; double f0 = 0.2;

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

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

// Initial conditions t[0] = a; x[0] = 1.0; y[0] = 0.0;

// Frequency range double w_ini = 0.1; double w_fin = 3.0; double w_pas = 0.01;

// Loop over different values for frequency for (double w = w_ini; w

// Amplitude double amplitude = sqrt(x[N - 1] * x[N - 1] + y[N - 1] * y[N - 1]);

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

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

return 0; }

How can I modify this code to include the x^3 dependency such that the Duffing equation becomes : and obtain the plot of amplitude versus w/w0