function f=rhsarc(t,y)
global xmin xmax ymin ymax
global fswitch arcswitch

global epsilon

n=length(y)/2;

y1=y(1:n);
y2=y(n+1:2*n);


% the following ODEs are the first three cases from Tutorial 9, Q1.
switch fswitch
  case 1
    f=[y1-y2;y1+y2-2*y1.*y2];
  case 2 
    f=[y2.*exp(y2);1-y1.^2];
  case 3
    f=[1-y1.*y2;(y1-1).*y2];
  case 4
%% add new cases here
  case 5
%% or here
  otherwise
%% or here!
end 

% work in arc-length
small=1e-4;
if ( arcswitch==1) 
   fmag=sqrt(f(1:n).^2+f(n+1:2*n).^2)+small; 
   f(1:n)=f(1:n)./fmag;
   f(n+1:2*n)=f(n+1:2*n)./fmag;
end
% stop equations blowing up outside plot area
pad=zeros(n,1);
y11=[y1;pad];
y22=[pad;y2];
index1=find(y11<xmin | y11>xmax);
index2=find(y22<ymin | y22>ymax);

f(index1)=0;
f(index1+n)=0;

f(index2)=0;
f(index2-n)=0;