MATLAB绘图：丁同仁ODE线素场绘制方法

clear all;
close all;
clc;
y = -5:0.05:5;%%xy的一维数据
x = -5:0.05:5;
f = @(x,y)  -(x^2+1)*(y^2-1)*x^(-1)*y^(-1);
[x,y] = meshgrid(x,y);%%生成[x,y]格点
%下面处理格点位置的向量，要使左半平面的向量相反
for i=1:1:201
    for j=1:1:201
        if x(i,j)<0
            vx(i,j)=-cos(atan(f(x(i,j),y(i,j))));
            vy(i,j)=-sin(atan(f(x(i,j),y(i,j))));
        else
            vx(i,j)=cos(atan(f(x(i,j),y(i,j))));
            vy(i,j)=sin(atan(f(x(i,j),y(i,j))));
        end
    end
end
quiver(x,y,vx,vy);%%生成向量图
hold on
%下面绘制方程某些特解的曲线
g1=@(x,y) y.^2-1+exp(-x.^2)./x.^2
fimplicit(g1,'r',[-5,5]);
hold on
g2=@(x,y) y.^2-1+0.5*exp(-x.^2)./x.^2
fimplicit(g2,'r',[-5,5]);
hold on
g3=@(x,y) y.^2-1+2*exp(-x.^2)./x.^2
fimplicit(g3,'r',[-5,5]);
hold on
g4=@(x,y) y.^2-1-1*exp(-x.^2)./x.^2
fimplicit(g4,'r',[-5,5]);
hold on
g5=@(x,y) y.^2-1-0.5*exp(-x.^2)./x.^2
fimplicit(g5,'r',[-5,5]);
hold on
g6=@(x,y) y.^2-1-2*exp(-x.^2)./x.^2
fimplicit(g6,'r',[-5,5]);
hold on
axis equal;