function integc(xmin,xmax,ymin,ymax,fnc)
% ----------------------------------------------------------------------
%     This is a set up routine for plotting integral curves.
%     It is interactive and calls ODE45 to get the solution
%     curves.
% ----------------------------------------------------------------------
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Written by Deborah Sulsky, Department of Mathematics and Statistics,
% University of New Mexico, Albuquerque, NM 87131
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Set plot parameters
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 values
        x0=0.0;
        xfinal=0.0;
        xf=2.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
             plot(v(1),v(2),'r*')
        end
% Plot forward or backward or both?
        opt=input('Enter: f for foward, b for backward, fb for both ','s');
%
if(opt=='f' | opt=='fb')
% Set initial data for forward integration
        x0=v(1);
        y0=v(2);
        xfinal=x0;
        nx=1;
        ny=1;
        nn=1;
% Solve the ODE in increments of xf until the boundary is reached
           while(nx==nn & ny==nn)
              xfinal=xfinal+xf;
              [x,y]=ode45(fnc,[x0,xfinal],y0,options);
              index_x=find(x>xmin & x<xmax);
              index_y=find(y>ymin & y<ymax);
              index=intersect(index_x,index_y);
              if(length(index)>0); 
                  plot(x(index),y(index),'r'); 
              end;
              nx=length(index_x);
              ny=length(index_y);
              nn=length(index);
           end;
%
end;
%
if(opt=='b' | opt=='fb')
%
% Set initial data for backward integration
        x0=v(1);
        y0=v(2);
        xfinal=x0;
        nx=1;
        ny=1;
        nn=1;
% Solve the ODE in increments of xf until the boundary is reached
           while(nx==nn & ny==nn)
              xfinal=xfinal-xf;
              [x,y]=ode45(fnc,[x0,xfinal],y0,options);
              index_x=find(x>xmin & x<xmax);
              index_y=find(y>ymin & y<ymax);
              index=intersect(index_x,index_y);
              if(length(index)>0);
                  plot(x(index),y(index),'r');
              end;
              nx=length(index_x);
              ny=length(index_y);
              nn=length(index);
           end;
end;
