T = 1;
NN = 2.^[2:8];
rk_error = zeros(size(NN));
ie_error = zeros(size(NN));
im_error = zeros(size(NN));
bd_error = zeros(size(NN));

function y = f(t,x)
  y = x - x^2;
end

for k = 1:length(NN)
  y = ode_rk4 ('f', [0,T], 1/2, NN(k));
  rk_error(k) = abs(y(NN(k)+1)-e/(1+e));
end

for k = 1:length(NN)
  y = ode_ieuler ('f', [0,T], 1/2, NN(k));
  ie_error(k) = abs(y(NN(k)+1)-e/(1+e));
end

for k = 1:length(NN)
  y = ode_bdf4 ('f', [0,T], 1/2, NN(k));
  bd_error(k) = abs(y(NN(k)+1)-e/(1+e));
end

for k = 1:length(NN)
  y = ode_imid ('f', [0,T], 1/2, NN(k));
  im_error(k) = abs(y(NN(k)+1)-e/(1+e));
end

loglog (NN, ie_error, "*;Implicit Euler;", ...
	NN, 10^-3 * NN.^-4, ';Slope 4;', ...
	NN, im_error, "x;Implicit Midpoint;", ...
	NN, 10^-2 * NN.^-2, ';Slope 2;', ...
        NN, rk_error, '*;RK4;', ...
	NN, 10^-2 * NN.^-1, ';Slope 1;', ...
	NN, bd_error, "o;BDF4;") 