// This function finds a range of values of x2 that bracket the
// change of sign for a value of x1, given the starting x2, the
// ending x2, the increment and the parameters for md_disc.sci.
 
function [r]=md_solve (x1, start, ends, increment, u1, u2, C1, C2)

   r = md_disc ([x1;start], u1, u2, C1, C2);
   vsign = sign (r);
   rold = r;

   for x2=start:increment:ends,
      r = md_disc ([x1;x2], u1, u2, C1, C2);
      if sign (r) ~= vsign then
         r=[x2-increment rold; x2 r];
         break;
      end;
      rold = r;
   end;