t0 = 0;
x0 = 1;
v0 = 0;
T = 40;
A = 0;
B = 0;
C = 0.3;
a = 1;
b = 1;
c = 0;
N = 10000;

epsilon = 1e-4;

f = @(t, x) [(a - c * x(2) ^ 2) * x(1) - b * x(2) - A * x(2)^ 3 + B * cos(C * t); x(1)];
f_b_plus = @(t, x) [(a - c * x(2) ^ 2) * x(1) - (b + epsilon) * x(2) - A * x(2)^ 3 + B * cos(C * t); x(1)];
f_b_minus = @(t, x) [(a - c * x(2) ^ 2) * x(1) - (b - epsilon) * x(2) - A * x(2)^ 3 + B * cos(C * t); x(1)];

s_value = sqrt(a*a - 4 * b);

x_correct_f = @(t) exp(t/2 * (a - s_value)) * (exp(t* s_value) * (s_value - a) + s_value + a) / (2 * s_value);
v_correct_f = @(t) -b * exp(t/2 * (a - s_value)) * ((exp(t* s_value) - 1));

t = linspace(t0, T, N + 1);
xvs = zeros(2, N + 1);
vs = zeros(1, N + 1);
xs = zeros(1, N + 1);
xvs(1,1) = v0;
xvs(2,1) = x0;

for i = 1:N
  xvs(:,i + 1) = schemat(f, t(i), xvs(:,i), t(i + 1), 1);
  vs(i + 1) = v_correct_f(t(i));
  xs(i + 1) = x_correct_f(t(i));
endfor


mean_abs_diff_x = max(abs(xvs(2, :) - xs))
mean_abs_diff_v = max(abs(xvs(1, :) - vs))
printf("mean abs diff with x %f,\nmean abs diff with v %f\n", mean_abs_diff_x, mean_abs_diff_v)