function [t,y] = get_points_IV(ode_fxns,num_points,Tspan,IV,h)
a = Tspan(1);
t = linspace(a,a+(num_points-1)*h,num_points);
y = zeros(4,1);
Phi = zeros(4,1);
Phi_star = zeros(4,1);
y(:,1) = IV;
for i = 1:num_points
    %----------------------------------------------------------------- 
    k_1 = ode_fxns(t(i),y(:,i));
    k_2 = ode_fxns(t(i)+2/7*h, y(:,i)+2/7*h.*k_1);
    k_3 = ode_fxns(t(i)+7/15*h,y(:,i)+77/900*h.*k_1+343/900*h.*k_2);
    k_4 = ode_fxns(t(i)+35/38*h,y(:,i)+805/1444*h.*k_1...
                   -77175/54872*h.*k_2+97125/54872*h.*k_3);
    Phi(:,i) = (79/490.*k_1+2175/3626.*k_3+2166/9065.*k_4);
    %-----------------------------------------------------------------
    k_5 = ode_fxns(t(i)+h,y(:,i)+h*Phi(:,i));
    Phi_star(:,i) = 229/1470*k_1+1125/1813*k_3+13718/81585*k_4+1/18*k_5;
    y(:,i+1) = y(:,i) + h*Phi_star(:,i);
   
end
    
end