function ivp

%  solve IVP using MATLAB routines
%         y' = f(t,y)  with   y(0) = y0        
%  where y = (y1, y2, y3 , ..., yn) is an n-vector

% clear all previous variables and plots
clear *
clf

% time interval
tmax=200;

%  initial values
y10=0.5; y20=0; 
y0=[y10 y20];

%  calculate solution using a MATLAB routine
%[t,y] = ode45(@rhs,[0 tmax],y0);  titleX='\lambda=1';
[t,y] = ode23s(@rhs,[0 tmax],y0);  titleX='ode23s';

%
%  phase plane plot
%
hold on
plot(y(:,1),y(:,2),'r')

% put a marker at the start position
plot(y(1,1),y(1,2),'.r','LineWidth',1,'MarkerSize',36)

%  draw an arrow indicating direction of motion
i=30;
ii=i+1;
arrowhead([y(i,1) y(ii,1)],[y(i,2) y(ii,2)],'r');

% commands to label each axes
xlabel('y_1-axis','FontSize',14,'FontWeight','bold')
ylabel('y_2-axis','FontSize',14,'FontWeight','bold')

% have MATLAB use certain plot options (all are optional)
box on
% add a title to plot
title(titleX,'FontSize',14,'FontWeight','bold')
% Set the fontsize to 14 for the plot
set(gca,'FontSize',14); 
% Set legend font to 14/bold                            		
set(findobj(gcf,'tag','legend'),'FontSize',14,'FontWeight','bold'); 

% get(gcf)
%set(gcf,'Position', [625 537 524 251]);

hold off

%  define f(t,y)
function dy=rhs(t,y)
dy=zeros(2,1);
lambda=1;
dy(1) = y(2);
%dy(2) = +y(1)-y(2);
dy(2) = lambda*(1-y(1)^2)*y(2)-y(1);
%dy(2) = -lambda*y(1)-y(1)^3-y(2);