function [t,y] = sRK5(ode_f,Tspan,IV,h)
    a = Tspan(1);
    b = Tspan(2);
    n = (b-a)/h;
    y = zeros(4,1);
    Phi = zeros(4,1);
    Phi_star = zeros(4,1);
    y(:,1) = IV;
    t = a:h:b;
    t_start=1;
    j = 0;
%     figure(2);
%     hold on
    for i = 1:n
    %----------------------------------------------------------------- 
    k_1 = ode_f(t(i),y(:,i));
    k_2 = ode_f(t(i)+2/7*h, y(:,i)+2/7*h.*k_1);
    k_3 = ode_f(t(i)+7/15*h,y(:,i)+77/900*h.*k_1+343/900*h.*k_2);
    k_4 = ode_f(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_f(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);
%     if (mod(t(i),11.124/4) == 0)   
%         j=j+1;
%         if (j == 1)
%             plot(y(3,(t_start:i)),y(4,(t_start:i)),'red')
%             legend('j=1');
%         elseif(j ==2)
%             plot(y(3,(t_start:i)),y(4,(t_start:i)),'blue');
%             legend('j=1','j=2')
%         elseif(j==3)
%             plot(y(3,(t_start:i)),y(4,(t_start:i)),'Color',[0 0.5 0]);
%             legend('j=1','j=2','j=3')
%         elseif(j==4)
%             plot(y(3,(t_start:i)),y(4,(t_start:i)),'magenta');
%             legend('j=1','j=2','j=3','j=4')
%         elseif(j==5)
%             plot(y(3,(t_start:i)),y(4,(t_start:i)),'cyan');
%             legend('j=1','j=2','j=3','j=4','j=5')
%         elseif(j==6)
%             plot(y(3,(t_start:i)),y(4,(t_start:i)),'Color',[.5 0 0]);
%             legend('j=1','j=2','j=3','j=4','j=5','j=6')
%         elseif(j==7)
%             plot(y(3,(t_start:i)),y(4,(t_start:i)),'black');
%             legend('j=1','j=2','j=3','j=4','j=5','j=6','j=7')
%         else
%             plot(y(3,(t_start:i)),y(4,(t_start:i)),'yellow');    
%             legend('j=1','j=2','j=3','j=4','j=5','j=6','j=7')
%         end
%         t_start = i+1;
%     end
     end
%     hold off
end