function twodim 
% integrates 2-d autonomous system, specifically the 
% near-Hamiltonian system considered in Examples 26 and 27
%
LL=2;
xmin=-LL;
xmax=LL;
ymin=-LL;
ymax=LL;
%
NN=80;
deltax=(xmax-xmin)/NN;
deltay=(ymax-ymin)/NN;
%
xx=[-NN/2:1:NN/2]*deltax;
yy=[-NN/2:1:NN/2]*deltay;
%
H=0.5*(yy'.*yy'*ones(size(xx))+ones(size(yy))'*(xx.*xx)) ...
    -(1/3)*ones(size(yy))'*(xx.*xx.*xx);
cintH=[H(1+NN/2,1) H(1+NN/2,1+NN/8) H(1+NN/2,1+NN/4) ...
    H(1+NN/2,1+3*NN/8) H(1+NN/2,1+NN/2) H(1+NN/2,1+5*NN/8) H(1+NN/2,1+3*NN/4) ...
    H(1+NN/2,1+7*NN/8) H(1+NN/2,1+NN)];
contour(xx,yy,H,cintH, ...
    'LineColor','k');
xlabel('x');
ylabel('y');
%hold on
%
epsilon=input('epsilon= ');
aa=input('a= ');
FunFcn='nearh';
t0=0;
tfinal=100;
pow = 1/3;
tol = 0.001;
%
t = t0;
hmax = (tfinal - t)/5;
hmin = (tfinal - t)/200000;
h = (tfinal - t)/1000;
%y0(1)=0.2;
%y0(2)=0.2;
[y0(1),y0(2)]=ginput(1);
y = y0(:);
Y=y*ones(1,2);
%
while (h>=hmin) && (abs(y(1))<xmax) && (abs(y(2))<ymax)
    if t + h > tfinal, h = tfinal - t; end
    % Compute the slopes
    s1 = feval(FunFcn, t, y, epsilon, aa);
    s2 = feval(FunFcn, t+h, y+h*s1, epsilon, aa);
    s3 = feval(FunFcn, t+h/2, y+h*(s1+s2)/4, epsilon, aa);

            % Estimate the error and the acceptable error
            delta = norm(h*(s1 - 2*s3 + s2)/3,'inf');
            tau = tol*max(norm(y,'inf'),1.0);

            % Update the solution only if the error is acceptable
            if delta <= tau
                t = t + h;
                y = y + h*(s1 + 4*s3 + s2)/6;

                % Update the plot
                Y = [y Y(:,1)];
                %line( ...
                %  'color','r', ...
                %  'Marker','.', ...
                %  'markersize',25, ...
                %  'erase','background', ...
                %  'xdata',y(1),'ydata',y(2));
                line(Y(1,:),Y(2,:),'Color','b','LineWidth',2);
                drawnow;
                pause(0.1)
            end

            % Update the step size
            if delta ~= 0.0
                h = min(hmax, 0.9*h*(tau/delta)^pow);
            end
   
end  %Main loop ...
end
%===============================================================================
function ydot = nearh(t,y,epsilon,aa)
%   specifies rhs of differential equation


%epsilon=0.4;
%aa=0.05;
A = [ 0 1 
     y(1)-1 epsilon*(aa-y(1))];
ydot = A*y;

end

