%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% You may want to declare a 2xM matrix to hold 
% your solution (so you can plot the phase plane)
% where M is the number of grid points (including
% the points t=0, and t=1).  The first row of Y
% hold the values of y1, and the second those of
% y2.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Y=zeros(2,M);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The following line assigns the initial values
% given in the problem to the first column in
% your solution matrix.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
Y(:,1)=[1 1]';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Modify your time advancement loop as follows.
% Where F is the function defined in the file F.m
% Doing it this way will compute the values of
% y1(i+1) and y2(i+1) simultaneously.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i=1:M-1

Y(:,i+1)=Y(:,i)+h*F(Y(:,i)+.5*h*F(Y(:,i),t(i)),t(i)+.5*h);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The following line should produce the phase
% plane.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

plot (Y(1,:),Y(2,:));