clear all; clc;

% Problem setting
t0 = 0;
y0 = 1;
T = 2;

f = @(t,y) t * y;

% Error for various stepsizes h
hs = 2.^(-(1:10));
for i = 1:10,
    h = hs(i);
    t = t0:h:T;
    y_true = exp(t.^2/2);
    y = nystroem(h, f, t0, y0, T);
    err(i) = norm(y_true - y, Inf);
end

% Plotting decay of error
figure(1), clf
loglog(hs, err,'-ro', hs, hs.^3, '--k','LineWidth',2)
xlabel('stepsize h','FontSize',14)
legend('error of nystroem','h^3')
set(gca,'FontSize',13)

% Regression for convergence order
polyfit(log(hs),log(err),1)
eoc(hs', err')
