errVSstep.m

function errVSstep(A, y0, t0, tf)
errors = [];
range = 1:20;
h = (tf-t0)./2.^range;

for i=range
    [approx, err] = eulerint(A,y0,t0,tf,2^i);
    errors = [errors, norm(err)];
end

loglog(h,errors);

end

eulerint.m

function [approx, err] = eulerint(A, y0, t0, tf, N)
approx=0;
h = (tf-t0)/N;
y_t = y0;
for i=0:N
    approx = approx + y_t;
    y_t = eulerstep(A,y_t,h);
end
approx = approx/N;

exact = (expm(A*tf)-expm(A*t0))/A*y0;
err = exact-approx;

eulerstep.m

function unew = eulerstep(A, y_t, h)
unew = y_t+h*A*y_t;