function twodim
% integrates 2-d autonomous system, specifically the 
% system considered in Example 20
%
LL=4;
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;
%
cla;
axis([xmin xmax ymin ymax]);
line([xmin xmax],[0 0],'Color','k');
line([0 0], [ymin ymax],'Color','k');
box on;
xlabel('x');
ylabel('y');
hold on
%
aa=input('a= ');
FunFcn='ff';
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, aa);
    s2 = feval(FunFcn, t+h, y+h*s1, aa);
    s3 = feval(FunFcn, t+h/2, y+h*(s1+s2)/4, 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 ...
hold off
end
%===============================================================================
function ydot = ff(t,y,aa)
%
% specifies rhs of differential equation in Example 20
%
A = [ 1+2*y(1)-y(1)*y(1)-y(2)*y(2) -1+aa*y(1) 
     1-aa*y(1) 1+2*y(1)-y(1)*y(1)-y(2)*y(2)];
ydot = A*y;
end

