%   This program plots the direction field or integral curves
%   for a first order system of equations
%                       dx/dt = f(x,y)
%                       dy/dt = g(x,y)
%
%   The user must supply the function defining the first order system
%
%   The syntax is
%            function yprime=f(t,z)
%   where z=[x,y] the point in the xy-plane and yprime is the corresponding
%   slope; i.e. yprime contains dx/dt and dy/dt at (x,y).
%   (yp(1)=f(x,y) and yp(2)=g(x,y))
%   This function must return a column vector for yprime (for ODE45)
%   t is time which must be passed along for ODE45 
%
%
%   The direction field is plotted at (nx)x(ny) equally spaced
%   points in the area defined by (xmin,xmax)x(ymin,ymax).
%
%   Parameters:
%         xmin,xmax   minimum and maximum values of x to use
%         ymin,ymax   minimum and maximum values of y to use
%         nx          number of points in the x-direction
%         ny          number of points in the y-direction
%         fncname     name of the function defining the system of 2 ODEs
%
%   Integral curves ask for the starting point and the direction of
%   integration.  Integration is carried out for 2 time units or until
%   the boundary is reached.  An option is given to continue the same
%   integral curve in the same direction or to reverse direction (or stop).
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Written by Deborah Sulsky, Department of Mathematics and Statistics,
%     University of New Mexico, Albuquerque, NM 87131
%
% Get the name of the function to use
  fncname=input('Enter the name of the file containing the function: ','s');
%
% Get the plotting window
    wind=input('Enter the plot size: [xmin,xmax,ymin,ymax] ');
    xmin=wind(1); xmax=wind(2); ymin=wind(3); ymax=wind(4);
%
% Set control for direction field
%      done=0  allows direction field to be plotted
%      afterward done=1 so the direction field can be plotted only once.
       done=0;
% Set control for execution to continue
       quitpp=0;
%
% Open figure windows
       h1=figure;
       h2=figure;
%
%
while(quitpp==0);
%
%
    % Determine action
    disp('Enter option:')
    disp('     d for plotting the direction field')
    disp('     i for plotting an integral curve')
    action=input('     q to quit \n','s');
%
    if(action=='d');
%
        %  If direction field already done, do not plot it again
        if (done==0)
%
        %  Get the number of points to use
        disp('The direction field is plotted at nx by ny points')
        N=input('enter the number of points: [nx,ny] ');
        nx=N(1); ny=N(2);
%
        %   Clear figure window
        clf
        %   Clear storage
        if(exist('x')) clear x; end;
        if(exist('y')) clear y; end;
        if(exist('X')) clear X; end;
        if(exist('Y')) clear Y; end;
        if(exist('F')) clear F; end;
%
        %   Determine the spacing between points
        hx=(xmax-xmin)/nx;
        hy=(ymax-ymin)/ny;
        %   Make a grid of x and y values, using this spacing
        x=xmin+hx/2:hx:xmax-hx/2;
        y=ymin+hy/2:hy:ymax-hy/2;
        [X,Y]=meshgrid(x,y);
%
        for i=1:length(x);
           for j=1:length(y);
              %Evaluate the slope of the direction field at (x(i),y(j))
              dummyt=1.0;
              yp=feval(fncname,dummyt,[x(i);y(j)]);
              %Scale the line segments so they all have the same length
              scale=sqrt(yp(1)*yp(1)+ yp(2)*yp(2));
              yp=yp/scale;
              % Save the values for plotting
              F(j,i,1)=yp(1);
              F(j,i,2)=yp(2);
           end
        end
%
        %   Plot the direction field using quiver
        figure(h1);
        quiver(X,Y,F(:,:,1),F(:,:,2),0.5);
        %   Set the plot frame
        axis([xmin,xmax,ymin,ymax]);
        %   Add axes
        hold on;
        plot([xmin,xmax],[0,0],'k');
        plot([0,0],[ymin,ymax],'k');
        %   Direction field done
        done=1;
%
        else;
%
        disp('Direction field already done!!')
        disp('')
        end;
%
%
%
    elseif(action=='i');
%
%
%
        % Integral Curves
        % Set plot parameters
        figure(h1);
        hold on;
        axis([xmin,xmax,ymin,ymax]);
        % Set options for solver
        oldopts=odeset;
        options=odeset(oldopts,'RelTol',1e-8,'AbsTol',1e-8,'NormControl','off');
        % Set default start values
        x0=0.0;
        xfinal=0.0;
        xf=2.0;
	cc=0;
%
        % Get initial value
        v=input('Enter initial value for the integral curve: [x0,y0] ');
        % Make sure v is in bounds
        if(v(1)< xmin | v(1)>xmax | v(2)< ymin | v(2)>ymax)
             error('initial data outside domain')
        else
             figure(h1);
             plot(v(1),v(2),'r*')
        end
        % Plot forward or backward or both?
		disp('For this integral curve:')
		disp('  Enter f for Forward')
        opt=input('  Enter b for Backward \n','s');
%
     while(cc==0);
        if(opt=='f')
           % Set initial data for forward integration
           x0=0.0;
           xfinal=x0;
           y0=v';
           % Solve the ODE in an increment of xf 
	      xf=abs(xf);
              xfinal=xfinal+xf;
              [x,y]=ode45(fncname,[x0,xfinal],y0,options);
              index_x=find(y(:,1)>xmin & y(:,1)<xmax);
              index_y=find(y(:,2)>ymin & y(:,2)<ymax);
              index=intersect(index_x,index_y);
              if(length(index)>0);
                  figure(h1);
                  plot(y(index,1),y(index,2),'r');
                  figure(h2)
                  plot(x(index),y(index,1),x(index),y(index,2))
                  legend('y_1(t)','y_2(t)')
                  hold on;
              end;
        end;
%
        if(opt=='b')
           % Set initial data for backward integration
           x0=0.0;
           xfinal=x0;
           y0=v';
           % Solve the ODE in an increment of xf 
	      xf=-abs(xf);
              xfinal=xfinal+xf;
              [x,y]=ode45(fncname,[x0,xfinal],y0,options);
              index_x=find(y(:,1)>xmin & y(:,1)<xmax);
              index_y=find(y(:,2)>ymin & y(:,2)<ymax);
              index=intersect(index_x,index_y);
              if(length(index)>0);
                  figure(h1);
                  plot(y(index,1),y(index,2),'r');
                  figure(h2)
                  plot(x(index),y(index,1),x(index),y(index,2))
                  legend('y_1(t)','y_2(t)')
                  hold on;
              end;
	  end;
%
          if(opt=='c')
          % Set initial data to continue integration
		  x0=xfinal;
          y0=y(length(y),:);
          % Solve the ODE in an increments of xf
              xfinal=xfinal+xf;
              [x,y]=ode45(fncname,[x0,xfinal],y0,options);
              index_x=find(y(:,1)>xmin & y(:,1)<xmax);
              index_y=find(y(:,2)>ymin & y(:,2)<ymax);
              index=intersect(index_x,index_y);
              if(length(index)>0);
                  figure(h1);
                  plot(y(index,1),y(index,2),'r');
                  figure(h2)
                  plot(x(index),y(index,1),x(index),y(index,2))
                  legend('y_1(t)','y_2(t)')
                  hold on;
              end;
	  end;
%
        if(opt=='r')
           % Set initial data for reverse integration
           x0=0.0;
           xfinal=x0;
           y0=v';
           % Solve the ODE in an increment of xf
              xf=-xf;
              xfinal=xfinal+xf;
              [x,y]=ode45(fncname,[x0,xfinal],y0,options);
              index_x=find(y(:,1)>xmin & y(:,1)<xmax);
              index_y=find(y(:,2)>ymin & y(:,2)<ymax);
              index=intersect(index_x,index_y);
              if(length(index)>0);
                  figure(h1);
                  plot(y(index,1),y(index,2),'r');
                  figure(h2)
                  plot(x(index),y(index,1),x(index),y(index,2))
                  legend('y_1(t)','y_2(t)')
                  hold on;
              end;
              end;
%
	  disp('Enter c to continue this integral curve')
	  disp('Enter r to reverse the direction')
	  opt=input('Enter s to stop \n','s');
	  if(opt=='s') cc=1; end;
     end;
%
elseif(action=='q')
     quitpp=1;
     close(h1);
     close(h2);
%
end;
end;
