ti = 0;
tf = 10*pi;

n = 1000;
h = (tf-ti)/(n-1); 

yi = [1,0];

yout = zeros(n,2);

t = ti;
for i = 1:n
    yout(i,:) = yi;
    yi = rk4(yi,t,h,@sho);
    t = t+h;
end