t0=0;  tN=20;    y0=[0;0.25];   h=0.001;

t = t0:h:tN;      N = length(t);    

j = 1;
for i = 1:N
      t1 = t0 + h;
      K1 = m_chap2_ex1_1_sub(t0, y0);
      K2 = m_chap2_ex1_1_sub(t0 + h/2, y0 + h*K1/2);
      K3 = m_chap2_ex1_1_sub(t0 + h/2, y0 + h*K2/2);
      K4 = m_chap2_ex1_1_sub(t0 + h, y0 + h*K3);
      y1 = y0 + (h/6)*(K1 + 2*K2 + 2*K3  + K4);
      yy1(j)=y1(1);  yy2(j)=y1(2);  
      t0=t1;   
      y0=y1;
      j=j+1;
end

t=0:h:tN;

subplot (121),   plot(t,yy1,t,yy2);  grid
subplot (122),   plot(yy1,yy2);      grid