#include <iostream>
#include <iomanip>

#include <math.h>

/****************************************************************/

// Since this is just an example and not a serious implementation
//   I'm just coding the matrix vector multiplication and the
//   matrix inverse folrmulas for 2x2 matrices explicitly.


inline void Ax_to_x(double A[2][2], double x[2])
{
   double x0 = A[0][0]*x[0]+A[0][1]*x[1];
   double x1 = A[1][0]*x[0]+A[1][1]*x[1];

   x[0]=x0; x[1]=x1;
}

double invertA(double A[2][2])
{
   double a00 = A[0][0]; double a01 = A[0][1];
   double a10 = A[1][0]; double a11 = A[1][1];

   double det = a00*a11 - a10*a01;

   if (det!=0)
   {
      A[0][0] = a11/det; A[0][1] = -a01/det;
      A[1][0] = -a10/det; A[1][1] = a00/det;      
   }

   return det;
}

/****************************************************************/

int main()
{
   double x[2] ={1,0};
   double A[2][2];
   A[0][0] = -1000; A[0][1] = 1;
   A[1][0] = 999;   A[1][1] = -2;

   double h = 0.01;

   // Implicit Euler Method propagation matrix (x_(n+1) = (I-hA)^{-1}x_n))
   double IE[2][2];
   IE[0][0] = 1-h*A[0][0];    IE[0][1] = 0-h*A[0][1];
   IE[1][0] = 0-h*A[1][0];    IE[1][1] = 1-h*A[1][1];
   invertA(IE);

   // Explicit Euler Method propagation matrix (x_(n+1) = (I+hA)x_n))
   double EE[2][2];
   EE[0][0] = 1+h*A[0][0];    EE[0][1] = 0+h*A[0][1];
   EE[1][0] = 0+h*A[1][0];    EE[1][1] = 1+h*A[1][1];
 
   for (int i=0; i<100; i++)
   {
      std::cout << i*h << " - " << "(" << x[0] << "," << x[1] << ") = ("
                << 0.999*exp(-1001*i*h)+0.001*exp(-i*h) << ","
                << -0.999*exp(-1001*i*h)+0.999*exp(-i*h) << ")" << std::endl;
      Ax_to_x(IE,x);
   }
}