function [x,y]=tdrknu(m,h,x,y,subl,subn)

  dx=h/m;
  for ii=1:m
    [x,y] = tdrknl(dx,x,y,subl,subn);
  end
