function [t,y,Phi,Phi_star,v] = Sys4ODEsRK4(ode_f,Tspan,IV,h,q)
  %output2 = fopen('output2','a');
  output3 = fopen('output3','a');
    a = Tspan(1);
    b = Tspan(2);
    n = (b-a)/h;
    five_terms=0;
    y = zeros(4,n);
    Phi = zeros(4,n);
    Phi_star = zeros(4,n);
    r = zeros(4,n);
    v = zeros(4,n);
    y(:,1) = IV;
    v(:,1) = 0;
    t = [a:h:b];
    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);
    %-----------------------------------------------------------------
    if(five_terms==0)
        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);
        r(:,i) = h^3*(Phi(:,i)-Phi_star(:,i));
        v(:,i+1) = v(:,i) + h*(Phi(:,i).*v(:,i)+r(:,i));
    else
        y(:,i+1) = y(:,i) + h*Phi(:,i);
    end
    if (t(i)==0|t(i)==0.2|t(i)==0.4|t(i)==0.6|t(i)==0.8)
    %    fprintf(output2,'\t %i \t %f \t %f \t %d \t %d \n',q,h,t(i),max(...
    %    [max(abs(y(3,i+1)-cos(pi/2*t(i+1)^2))) max(abs(y(4,i+1)-sin(pi/2*t(i+1)^2)))])...
    %    ,h*max(Phi(:,i)-Phi_star(:,i)));
    fprintf(output3,'\t %i \t %f \t %f \t %d \t %d \n',q,h,t(i),max(...
        [max(abs(y(3,i+1)-cos(pi/2*t(i+1)^2))) max(abs(y(4,i+1)-sin(pi/2*t(i+1)^2)))])...
        ,h^3*max(v(:,i+1)));
        if(t(i) == 0.8) 
            fprintf(output3,'--------------------------------------------------------------\n');
        end
    end
%    end
    
end