function [UU,VV]=dirfieldsolns(ff,x,y,y0,h);
% DIRFIELDSOLNS  plots a direction field, and superimposes solutions
%                which start at left edge, computed by Euler method
%
%  Example:  Assume you want the direction field of
%    dy
%    --  = x + 2*sin(y)      for -1<=x<=1, -3<=y<=3
%    dx
%  and you want solutions with initial values y(-1) = 0, 0.5, 2.
%
%  Write a function file "gg.m" as follows:
%     function dy = gg(x,y);
%     dy= x + 2*sin(y);
%
%  Then run:
%  >> dirfieldsolns('gg',-1:.2:1,-3:.5:3,[0 0.5 2]',0.01);

[XX,YY]=meshgrid(x,y);
ANG=atan(feval(ff,XX,YY));

UU=cos(ANG);
VV=sin(ANG);

quiver(XX,YY,UU,VV,.25,'k.');
hold on;
quiver(XX,YY,-UU,-VV,.25,'k.');

xeuler = min(x):h:max(x);
for n = 1:length(y0)  % run through initial values
  yeuler = zeros(size(xeuler));
  yeuler(1) = y0(n);
  for k = 2:length(xeuler)  % step from left to right using h
    yeuler(k) = yeuler(k-1) + h * feval(ff,xeuler(k-1),yeuler(k-1));
  end
  plot(xeuler,yeuler)
end
hold off