% styveul, styvt system   
%          du1/dt=a1*u1+a2*u2,  du2/dt=(a1+1)*u1+(a2-1)*u2
  clear, clf
  ustart=[-2 2]';  tslut=1;
  a1=-200; a2=120;  
  A=[a1  a2;  a1+1  a2-1];
% Kolla egenvärdena med eig(A): 
  lambda=eig(A)                 % resultat:  -80 och -1
  disp('Eulers metod:')
% stabil om dt < 2/max(abs(lambda)) alltså  dt<2/80 = 0.025
  tidssteg=[0.05 0.025 0.02 0.01];
  disp('     dt        u1        u2 ')
  for nr=1:4 
    dt=tidssteg(nr); u=ustart;  t=0;  U=u';  T=t;
    while t<tslut-dt/2
      u=u+dt*A*u; t=t+dt;  
      U=[U; u'];  T=[T; t];
    end
    disp([dt u'])
    Umax=max(U); Umin=min(U);
    subplot(3,2,nr), plot(T,U), title(['Euler, dt= ' num2str(dt)])
    axis([0 1 min(Umin) max(Umax)]),  drawnow
  end
  disp('Implicit Euler:')
  tidssteg=[0.1 0.05];
  for nr=1:2
    dt=tidssteg(nr); u=ustart;  
    t=0; U=u'; T=t;
    C=eye(2)-dt*A;                   % (I-dt*A)u(t+dt)=u(t)
    while t<tslut-dt/2
      u=C\u;     t=t+dt;  
      U=[U; u']; T=[T; t];
    end
    disp([dt u'])
    subplot(3,2,4+nr), plot(T,U), title(['implicEuler, dt= ' num2str(dt)])
    axis([0 1 -2 10]),  drawnow
  end



