function hfig = vector_field(f,y0,t)
%vector_field Calculate and plot the vector field of f
%
%   INPUT:  f       Function handle of t and y, calculating the derivatives
%                   dy/dt
%           y0      Vector of different initial conditions to solve for
%           t       Vector of time points where the vector field is plotted
%
%   OUTPUT: hfig    Handle to the figure that is created



Y = zeros(length(y0),length(t));

for i = 1:length(y0)
    [T,Y(i,:)] = ode45(f,t,y0(i));
end

T = repmat(t,length(y0),1);

dy = f(T,Y);
U = ones(size(Y));
V = dy;
L = sqrt(U.^2+V.^2);

hfig = quiver(T,Y,U./L,V./L);
axis tight