% plot direction field for first order ODE y' = f(t,y)
function dirfield
              
f = @(t,y) -y^2+t % f = @(independent_var, dependent_var)
    % mathematical expression y in terms of t

figure;             
dirplotter(f, -1:.2:3, -2:.2:2) % dirfield(f, t_min:t_spacing:t_max, y_min:y_spacing:y_min)
    % using t-values from t_min to t_max with spacing of t_spacing
    % using y-values from y_min to t_max with spacing of y_spacing
end

function dirplotter(f,tval,yval) %  plots direction field using aboved defined (f,tval,yval)
[tm,ym]=meshgrid(tval,yval); % creates a grid of points stored in matrices in the (t,y)-plane
  dt = tval(2) - tval(1); % determines t_spacing and y_spacing based on tval
  dy = yval(2) - yval(1); % determines t_spacing and y_spacing based on yval
  fv = vectorize(f); % prevents matrix algebra; if f is a function handle fv is a string
  if isa(f,'function_handle') % if f is function handle
  fv = eval(fv); % turn fv into a function handle from a string 
  end
  yp=feval(fv,tm,ym); % evaluates differential equation fv from tm and ym from meshgrid
  s = 1./max(1/dt,abs(yp)./dy)*0.35;
  h = ishold; % define h as ishold which returns 1 if hold is on, and 0 if it is off
  quiver(tval,yval,s,s.*yp,0,'.r'); hold on; % plots vectors at each point w direction based on righthand side of differential equation
  quiver(tval,yval,-s,-s.*yp,0,'.r');
  if h % h above defined as ishold
  hold on % retain current plot 
  else % refers to all not h 
  hold off % resets axes properties to default
  end
  axis([tval(1)-dt/2,tval(end)+dt/2,yval(1)-dy/2,yval(end)+dy/2]) % sets the scaling for axes on current graph to ends of tval and yval   
end    

% Simple direction field plotter
clear;
[T Y]=meshgrid(-2:0.2:2,-2:0.2:2);
dY=cos(3*T)+Y./T;
dT=ones(size(dY));
scale=sqrt(1+dY.^2);
figure;
quiver(T, Y, dT./scale, dY./scale,.5,'Color',[1 0 0]);
xlabel('t');
ylabel('y');
title('Direction Field for dy/dt');
axis tight;
hold on;

% change color
% do multiple examples, x.^2-y.^2; use T, etc

% Simple direction field plotter
clear;
[T Y]=meshgrid(-2:0.2:2,-2:0.2:2);
dY=T.^2-Y.^2;
dT=ones(size(dY));
scale=sqrt(1+dY.^2);
figure;
quiver(T, Y, dT./scale, dY./scale,.5,'Color',[0 1 0]);
xlabel('t');
ylabel('y');
title('Direction Field for dy/dt');
axis tight;
hold on;

			Complete