function [t,Y]=ImpEulerVec(F,interval, y0, m, const)
%function [t,Y]=EulerVec(F, interval, y0, m, const)
%  Solves the IVP (Improved Euler):  Y'=F(t,Y), y(interval(1))=y0 using m panels (m+1)
%  points.  Returns a matrix Y so that y(t) is Y(:,t) (that is, y is a
%  column vector).
%
%  The extra argument, const, is for parameters that need to be fed into
%  the ODE, as Y'=F(t,y,const)

t=linspace(interval(1),interval(2),m+1);
h=t(2)-t(1);
dim=length(y0);

Y=zeros(dim,m+1);
Y(:,1)=y0(:);  %This notation forces y0 to be a column vector

for j=1:m
    M1=F(t(j),Y(:,j),const);
    K1=Y(:,j)+h*M1;
    M2=F(t(j+1),K1,const);
    Y(:,j+1)=Y(:,j)+(h/2)*(M1+M2);
end