% initial values
t0 = 0;
y0 = [3; 1];

% rhs of ode
f = @(t,y) [-y(1)*y(2)^2*sin(pi*t); 5*y(2)*sin(pi*t)];

% setup butcher tableau for Cash-Karp method (5th order)
A = [[0,0,0,0,0,0];[1/5,0,0,0,0,0];[3/40,9/40,0,0,0,0];[3/10,-9/10,6/5,0,0,0];[-11/54,5/2,-70/27,35/27,0,0];[1631/55296,175/512,575/13824,44275/110592,253/4096,0]];
b = [2825/27648, 0, 18575/48384, 13525/55296, 277/14336, 1/4];
c = [0, 1/5, 3/10, 3/5, 1, 7/8];

% compute ivp solution
h = 1e-2;
t_limit = 2;
[times, steps] = run_rk(A, b, c, t0, y0, f, t_limit, h);
disp('y(t_limit) = ')
disp(steps(:,end))