clear
f = inline('2 * y * t','y','t')


Ydv =[0 0 0 0 ];
Yr = [];
Tr = [];

	t0 = 0.1;
    y0 = 1;
	h = 0.05;
	n = 20;
    t1 = 0;
    
     if (n <=0)
        n = abs(t1 - t0)/h;
    end
	
        
        y = y0;
	    t = t0;


	

	for j = 1:4
	k1 = feval(f,y,t);
	k2 = feval(f,y+   (h*k1*0.5) ,   t+(h*0.5));
        k3 = feval(f,y+   (h*k2*0.5) ,t+(h*0.5)); 
        k4 = feval(f,y+   (h*k3) ,       t+h);
        y = y + (h/6)*(k1 + 2*k2 + 2*k3 + k4);
	t = t+h;
	Ydv(j) = feval(f,y,t);
	Yr(j) = y;
	Tr(j) = t;
        end
	

	for j = 5:n
         u = y + (h/24)*(55*Ydv(abs(rem(j-2,4))+1) - 59*Ydv(abs(rem(j-3,4))+1) + 37*Ydv(abs(rem(j-4,4))+1) - 9*Ydv(abs(rem(j-5,4))+1));
     	t = t+h;
        Ydv(abs(rem(j-1,4))+1) = feval(f,u,t);
        y = y + (h/24)*(9*Ydv(abs(rem(j-1,4))+1) + 19*Ydv(abs(rem(j-2,4))+1) - 5*Ydv(abs(rem(j-3,4))+1) + Ydv(abs(rem(j-4,4))+1));
	Yr(j) = y;
	Tr(j) = t;        
end
	
R = [Tr ; Yr];
R'
       
	
	

