% Direction field for autonomous 2 by 2 systems m-file (script) % 3.2.1998 by HA (modified from Bäckström's 1 eq. case % Requires the system m-file to be written to work for an initial value % matrix (with 2 rows) instead of a column vector. This form is perfectly % good for the ode-solvers % clear variables; hold off; xmax=10; ymax=3; delx=xmax/20; dely=ymax/20; % *** problem dependent x=-4:delx:xmax; y=-ymax:dely:ymax; [X,Y]=meshgrid(x,y); Z=sys5(0,[X(:),Y(:)]'); DX=Z(1,:);DY=Z(2,:); DX=reshape(DX,size(X));DY=reshape(DY,size(Y)); L=sqrt(DX.^2+DY.^2); % Initial lengths L=L+1e-10; DX1=DX./L; DY1=DY./L; % Unit length for all arrows figure(1); quiver(X,Y,DX1,DY1,'g'); grid on;title('Direction field'); % axis([x1 x2 y1 y2]);
After having played with meshgrid, reshape and Mat(:) you should now be able to read the above dir2-code as an open book.
function u=sys5(t,y) % %global k c; k=1;c=1/2; % Make these local while testing u(1,:)=y(2,:); u(2,:)=-c*y(2,:)-k*sin(y(1,:));
hold on init=ginput(1);[T,Y]=ode45('sys5',[0,10],init);plot(Y(:,1),Y(:,2),'g') axis([5 7 -1 1]) init=ginput(1);[T,Y]=ode45('sys5',[0,20],init);plot(Y(:,1),Y(:,2),'r') for i=1:5; init=ginput(1);plot(init(1),init(2),'o');[T,Y]=ode45('sys5',[0,20],init); plot(Y(:,1),Y(:,2),'r') end;