f = inline('y*(y - 1)','y','t')
   
    t0 = 0;
    y0 = 0.5;
	h = 0.5;
	n = 6;
    t1 = 0;
    
    if (n <=0)
        n = abs(t1 - t0)/h;
    end
    
    Yr =[];
    Tr =[];
    
     y = y0;
	 t = t0;
     
     for j = 1:n
	 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;
	 Yr(j) = y;
	 Tr(j) = t;
	end
        
	R = [Tr ; Yr];
    R'