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); }