function [t,y]=RK(fun,tspan,y0,n,A,b,c)

        h=(tspan(2)-tspan(1))/n;
         y(:,1)=y0;
         t(1)=tspan(1);
         s=length(b);
         for j=1:n
 $              k(:,1)=feval(fun,t(j),y(:,j));
                        for i=2:s
                                k(:,i)=feval(fun,t(j)+c(i)*h,y(:,j)+h*k(:,1:i-1)*A(i,1:i-1)');
               end
         end
     2. function dy=pred_prey(t,y)
         function s=cos(x)
         k=1;
         a=2/3;
         d=4/3;
         r(t)=int(s(x^2),x,0,t);
         mu(t)=13/20-(3/5)*exp(-(3/x));
         dy(1)=(y(1)+k)*r(t)-a*y(1)*y(2);
         dy(2)=-mu(t)*y(2)+d*y(1)*y(2);
         dy=dy';
     3. function []=driver_pred_prey()
        close all
       tspan=[0 100];
       n=1000;
       fun='pred_prey';
       y0=[0 -1]';
       A=[0 0 0 0; 1/2 0 0 0 ; 0 1/2 0 0; 0 0 1 0];
       b=[1/6 1/3 1/3 1/6]';
       c=[0 1/2 1/2 1]';
       $$[t,y]=RK('pred_prey', tspan, y0, n, A, b, c);