clear s t x0 y0 a bsyms s t%s->x ,t对应y%f=sin(s)*sin(t);%f=t-s^2;f=-2*s*t/(1+s^2);% f=cos(s);a=32.0;b=32.0;x0=-16;y0=-16;m=40;n=40;h1=a/m;h2=b/n;hold %下面的2/3是线段缩放因子,需要在x,y两个方向同时乘上for I=1:m   s=x0+(I-1)*h1;   for j=1:n      t=y0+(j-1)*h2;      d=eval(f); %d表示斜率,如果解曲线族经过点(s,t)那么这点的斜率是d      y1=t+2/3*h1*d; %以x的增量计算y      if abs(y1-t)>2/3*h2 %(y1-t)=t+2/3*h1*d -t=2/3 * h1 * d         x1=s+1/d*h2*2/3; %以y的增量计算x,所以要乘以1/d;         plot([s,x1],[t,t+h2*2/3])      else          plot([s,s+h1*2/3],[t,y1])          %plot([1-3],[1,2])画一条直线，x为1到3,y为1到2          %在s(相当于x)的点加上h1(小步增量*2/3)          %注意上面的y1中的h1也乘于2/3,以缩短斜线长度      end   endend