clear all;
y0 = 1;
x = linspace(0, 2*pi, 300);
ysol = lsode("f", y0, x);
% for some reason the solution is put in a horizontal array instead of vertical. Transpose it:
y = ysol';
hold off;
plot(x, y);
 
% pretty plot axis scales:
axis tight
ylim auto
xlabel('Time x');
ylabel('Solution y(x)');

%% b:
hold on
yanal = sin(x).^2 - cos(x);
plot(x, yanal, 'r-');
ynum = diff(y) ./ diff(x);
plot(x(1:299)+(x(2)-x(1))/2, ynum, 'g-');
% (x(2)-x(1))/2 is to shift the curve to the right slightly to put derivative between two points