v0 = [ 0; 0; 10; 10 ];
ax = 0;
ay = -9.81;
f = @(v) [ v(3); v(4); ax; ay ];

len = 10;

result = zeros(size(v0,1), len+1);
h = zeros(1, len-1);
e = zeros(1, len-1);

for i = 0:len
	h(i+1) = 2^-i
	V = euler(f, v0, h(i+1), T);
	result(:,i+1) = V(:,end);
end

% calculation of error:
for i = 0:len
	e(i+1) = norm(result(:,i+1)-result(:,len+1));
end

loglog(h, e)