clear
g = '2*t*y';
	function out = f(y,t)
	 out = evstr(g)
	endfunction


Ydv =[0 0 0 0 ]
Yr =[];
Tr =[];
	t0 = 0.1;
        y0 = 1;
	h = 0.05;
	n = 20;
	
        y = y0;
	t = t0;


	

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

	for i = 5:n
        u = y + (h/24)*(55*Ydv(pmodulo(i-2,4)+1) - 59*Ydv(pmodulo(i-3,4)+1) + 37*Ydv(pmodulo(i-4,4)+1) - 9*Ydv(pmodulo(i-5,4)+1));

     	t = t+h;
        Ydv(pmodulo(i-1,4)+1) = f(u,t);
        y = y + (h/24)*(9*Ydv(pmodulo(i-1,4)+1) + 19*Ydv(pmodulo(i-2,4)+1) - 5*Ydv(pmodulo(i-3,4)+1) + Ydv(pmodulo(i-4,4)+1));
 	
	Yr(i) = y;
	Tr(i) = t;

	end

        R = [Tr Yr]