Xs = 5;
Xs_dot = 0;

%

m = 100;
k = 1000;
a = 400;
B = 200;

%

x0 = [-5; 0];

%

dxdt = @(t, x) [x(2); (B * (Xs_dot - x(2)) + k * (Xs - x(1)) + a * (Xs - x(1))^3) / m]

%

[t, x] = ode45(dxdt, [0, 5], x0)

%

plot(t, x(:, 1))

%

xk = linspace(0, 5);
fk = k * xk + a * xk.^3;
plot(xk, fk)