function test_phase_diagram
a = 1; b = 4; dx = 1/4;
c = -3; d = 3; dy = 1/6;
[x,y] = meshgrid(a:dx:b,c:dy:d);                                
r = (x.^2+y.^2).^(1/2);
v = sin(r); w = cos(r);
quiver(x,y,v,w,'c'); hold on;
v0 = 1.5; w0 = 2; 
streamline(x,y,v,w,v0,w0); hold off;