import graph_settings;
real stepfraction=0.05;

picture slopefield(real f(real,real), pair a, pair b,
                  int nx=nmesh, int ny=nx,
                  real tickfactor=0.5, pen p=currentpen)
{
 picture pic;
 real dx=(b.x-a.x)/nx;
 real dy=(b.y-a.y)/ny;
 real step=0.5*tickfactor*min(dx,dy);

 for(int i=0; i <= nx; ++i) {
   real x=a.x+i*dx;
   for(int j=0; j <= ny; ++j) {
     pair cp=(x,a.y+j*dy);
     real slope=f(cp.x,cp.y);
     real mp=step/sqrt(1+slope^2);
     draw(pic,(cp.x-mp,cp.y-mp*slope)--(cp.x+mp,cp.y+mp*slope),p);
   }
 }
 return pic;
}

picture slopefield(real f(real), pair a, pair b,
                  int nx=nmesh, int ny=nx, pen p=currentpen)
{
 return slopefield(new real(real x, real y) {return f(x);},a,b,nx,ny,p);
}

path curve(pair c, real f(real,real), pair a, pair b)
{
 real step=stepfraction*(b.x-a.x);
 real halfstep=0.5*step;
 real sixthstep=step/6;

 path follow(real sign) {
   pair cp=c;
   guide g=cp;
   real dx,dy;
   real factor=1;
   do {
     real slope;
     pair S(pair z) {
       slope=f(z.x,z.y);
       return factor*sign/sqrt(1+slope^2)*(1,slope);
     }
     pair S3;
     pair advance() {
       pair S0=S(cp);
       pair S1=S(cp+halfstep*S0);
       pair S2=S(cp+halfstep*S1);
       S3=S(cp+step*S2);
       pair cp0=cp+sixthstep*(S0+2S1+2S2+S3);
       dx=min(cp0.x-a.x,b.x-cp0.x);
       dy=min(cp0.y-a.y,b.y-cp0.y);
       return cp0;
     }
     pair cp0=advance();
     if(dx < 0) {
       factor=(step+dx)/step;
       cp0=advance();
       g=g..{S3}cp0{S3};
       break;
     }
     if(dy < 0) {
       factor=(step+dy)/step;
       cp0=advance();
       g=g..{S3}cp0{S3};
       break;
     }
     cp=cp0;
     g=g..{S3}cp{S3};
   } while (dx > 0 && dy > 0);
   return g;
 }

 return reverse(follow(-1))&follow(1);
}

path curve(pair c, real f(real), pair a, pair b)
{
 return curve(c,new real(real x, real y){return f(x);},a,b);
}