%-------------------------------------------------------
function [x,y]=tdrknl(dx,x,y,tdfunl,tdfunn)
  w1=feval(tdfunn,x,y);
  w2=y+(dx/2)*w1;
  w3=y+(dx/6)*w1;
 
  y =feval(tdfunl,x,dx/2,y);
  w2=feval(tdfunl,x,dx/2,w2);
  w3=feval(tdfunl,x,dx/2,w3);  
  w1=feval(tdfunn,x+dx/2,w2);

  w2=y+(dx/2)*w1;
  w3=w3+(dx/3)*w1;

  w1=feval(tdfunn,x+dx/2,w2);

  w2=y+dx*w1;
  w3=w3+(dx/3)*w1;

  w2=feval(tdfunl,x+dx/2,dx/2,w2);
  w3=feval(tdfunl,x+dx/2,dx/2,w3);  
  w1=feval(tdfunn,x+dx,w2);  

  x=x+dx;
  y=w3+(dx/6)*w1;
