斜率场--微分方程图形解

clear s t x0 y0 a b
syms 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
   end
end