
function rf(fun,x0,x1,y0,y1)
% rf skapar matriser X och Y som kan användas för att rita riktningsfältet
% till en differentialekvation y' = f(x,y) och plottar sedan
% riktningsfältet.
% Funktionen f(x,y) skall vara definierad i en funktionsfil eller med inline.
% [X,Y]=riktningsfalt('funktion',x0,x1,y0,y1) ger (2,420)-matriser X och Y.
% Vektorerna erhålls som de riktade sträckorna från (X(1,k),Y(1,k)) till  
% (X(2,k),Y(2,k)). Riktningen är beräknad i mittpunkten på sträckan.

N=20;
h=(x1-x0)/(N+1); k=(y1-y0)/N; [x,y]=meshgrid(x0:h:x1,y0:k:y1); 
% här genereras två (20,21)-matriser x och y som ger koordinaterna
% i ett rutnät med 420 punkter.

x=x(:)';y=y(:)';
% x(:)' och y(:)' är radmatriser med de 420 elementen i x resp y
% uppräknade kolonnvis.

yprim=feval(fun,x,y); 
% feval(fun,x,y) är en (20,21)-matris vars element är värdet av y' i rutnätet.

theta=atan(yprim);
% atan(feval(fun,x,y)) är tangentens riktningsvinkel.
 
r=0.4*min(h,k);u=r*cos(theta);v=r*sin(theta);
% (u(i),v(i)) är vektor med denna riktning och belopp=r.

X=[x-u;x+u]; Y=[y-v;y+v];
plot(X,Y)
