%---------------------------------------------------------  
% ルンゲ−クッタ法による時間発展(下位ルーチン)
%---------------------------------------------------------
function [x,y]=rk4l(dx,x,y,f)
  w1=feval(f,x,y);
  w2=y+(dx/2)*w1;
  w3=y+(dx/6)*w1;
  w1=feval(f,x+dx/2,w2);
  w2=y+(dx/2)*w1;
  w3=w3+(dx/3)*w1;
  w1=feval(f,x+dx/2,w2);
  w2=y+dx*w1;
  w3=w3+(dx/3)*w1;
  w1=feval(f,x+dx,w2);  
  x=x+dx;
  y=w3+(dx/6)*w1;