A = [-501 -499; -499 -501];
 x0 = [1;2];

 fprintf('matrix A:\n')
 disp(A)
 input('hit return')
 
 fprintf('\neigenvalues of A: \n')
 disp(eig(A))
 input('hit return')

 fprintf('\ninitial condition: \n')
 disp(x0)
 
 input('hit return')

% compute and plot the exact solution
 t = linspace(0,.025,200);
 xact = zeros(2,length(t));
 for k=1:length(t)
   xact(:,k) = expm(t(k)*A)*x0;
 end
 figure(1), clf
 set(gca,'fontsize',14)
 plot(t,xact(1,:),'b-','linewidth',2), hold on
 plot(t,xact(2,:),'r:','linewidth',2)
 legend('x_1(t)','x_2(t)')
 xlabel('t')
 ylabel('x(t)')
 title('Exact solution')
 input('hit return')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute the exact solution
 tend = .5;
 t = linspace(0,tend,200);
 xact = zeros(2,length(t));
 for k=1:length(t)
   xact(:,k) = expm(t(k)*A)*x0;
 end

 while 1
% solve with forward/backward Euler 
    method = input('\n Forward Euler (1) or Backward Euler (2) ? ');
    dt = input    ('\n                                 Delta t ? ');

    figure(1), clf
    set(gca,'fontsize',14)
    semilogy(t,abs(xact(1,:)),'b-','linewidth',2), hold on
    semilogy(t,abs(xact(2,:)),'r:','linewidth',2)
    xlabel('t')
    ylabel('|x(t)|')

    xk = zeros(2,ceil(tend/dt)+1);
    xk(:,1) = x0;
    if method==1,
        for k=1:ceil(tend/dt)
            xk(:,k+1) = xk(:,k) + dt*A*xk(:,k);     % forward Euler
        end
        title('Forward Euler')
    else
        for k=1:ceil(tend/dt)
            xk(:,k+1) = (eye(2)-dt*A)\xk(:,k);      % backward Euler
        end
        title('Backward Euler')
    end
    semilogy(dt*[0:ceil(tend/dt)],abs(xk(1,:)),'c.-','linewidth',2,'markersize',25)
    semilogy(dt*[0:ceil(tend/dt)],abs(xk(2,:)),'m.:','linewidth',2,'markersize',25)
    legend('xact_1','xact_2','eul_1','eul_2','Location','Best')
    drawnow
 end