g = 1;
f = @(t,x,eps) (eps*x-g*x.^3);
fx = @(t,x,eps) eps-3*g*x.^2;
M = 30;
dt = 0.001;
t = 0:dt:10;
x_lim = zeros(1,M);
eps = linspace(0,10,M);
for l=1:M
x = zeros(size(t));
x(1) = 0.1;
for i = 1:length(t)-1
x(i+1) = x(i)+dt*f(t(i),x(i),eps(l))+0.5*dt^2*(fx(t(i),x(i),eps(l)).*f(t(i),x(i),eps(l)));
end
x_lim(l) = x(end);
end
plot(eps,x_lim)
fitlm(log(eps(2:end)),log(x_lim(2:end)))