#include <math.h>
#include <fstream.h>
#include <iomanip.h>
#include <vector>


template <class Type>
void tri( int matrix_size, std::vector<Type> A, std::vector<Type> D, std::vector<Type> C, std::vector<Type> B, std::vector<Type> X)
{
 int i;
 Type xmult;

  for (i = 2; i <= matrix_size; i++)
  {
    xmult = A[i-1]/D[i-1];
    D[i] -= xmult * C[i-1];
    B[i] -= xmult * B[i-1];
  }
  X[matrix_size] = B[matrix_size]/D[matrix_size];
  for (i = matrix_size -1; i>=1; i--)
   X[i] = (B[i]- C[i]* X[i+1]) / D[i];

} 


template <class Type>
void seidel(Type ax, Type ay, int nx, int ny, Type h, int itmax, std::vector< std::vector<Type> >& u)
{
  int i, j, k;
  Type x, y, v;

  for (k = 1; k <= itmax; k++)
    {
      for (j = 1; j <= ny - 1; j++)
	{
	  y = ay + j*h;
	  for (i = 1; i <= nx - 1; i++)
	    {
	      x = ax + i*h;
	      v = u[i+1][j] + u[i-1][j] + u[i][j+1] + u[i][j-1];
	      u[i][j] = (v - h*h*g3(x,y)) / (4 - h*h*f3(x,y));
	    }
	}
    }
}