#include <iostream>
#include <iomanip>
#include <fstream>
#include <stdlib.h>

using namespace std;

const int A = 3;
const int N1 = 20;
const int N2 = 20;
const double phi = 0;
const double Eps = 0.001;
const double l1 = 1;
const double l2 = 1;
const double h1 = l1 / N1;
const double h2 = l2 / N2;

double k1(double x1, double x2) {
    return (x1 * x1 + x2 * x2 + 10);
}
double k2(double x1, double x2) {
    return (5*(x1 * x1 + x2 * x2)  + 20);
}
double p1(double x1, double x2) {
    return 0/*(-2*(x1 + 5*x2))*/;
}
double p2(double x1, double x2) {
    return 0/*(-1*(x1 * x1 + 5*x2*x2))*/;
}
double q(double x1, double x2) {
    return 0 /*(5*x1*x1 + x2*x2 + 30)*/;
}
double func(double x1, double x2){
    return(A * exp(x1 + x2));
}

double Af(double *mas[], double *x[] , int i , int j) {
    double T = (k1(x[0][i] + h1/2, x[1][j]) * (mas[i+1][j] - mas[i][j])) / (h1*h1)-
        (k1(x[0][i] - h1/2, x[1][j]) * (mas[i][j] - mas[i-1][j])) / (h1*h1)+
        (k2(x[0][i], x[1][j] + h2/2) * (mas[i][j+1] - mas[i][j])) / (h2*h2)-
        (k2(x[0][i], x[1][j] - h2/2) * (mas[i][j] - mas[i][j-1])) / (h2*h2)+
        (p1(x[0][i], x[1][j]) * (mas[i+1][j] - mas[i-1][j])) / (2*h1)+
        (p2(x[0][i], x[1][j]) * (mas[i][j+1] - mas[i][j-1])) / (2*h2)-
        q(x[0][i], x[1][j]) * mas[i][j];
    cout << T<< endl;
    system("pause");
    return T;
}
int main() {
    setlocale(LC_ALL, "");
    double** x = new double*[2];
    x[0] = new double[N1 + 1];
    x[1] = new double[N2 + 1];
    for (int i(0); i <= N1; i++) {
        x[0][i] = i * h1;
    }
    for (int j(0); j <= N2; j++) {
        x[1][j] = j * h2;
    }
    double** y = new double*[N1+1];
    double** u = new double*[N1+1];
    for(int i(0); i<=N1; i++) {
        y[i] = new double[N2+1];
        u[i] = new double[N2+1];
        for(int j(0); j<=N2; j++) {
            y[i][j] = 0;
            u[i][j] = 0;
        }
    }

    for(int i(0); i<=N1; i++) {
        for(int j(1); j<N2; j++) {
            y[i][0] = A * exp(x[0][i]);
            y[i][N2] = A * exp(l2+x[0][i]);
            y[0][j] = A * exp(x[1][j]);
            y[N1][j] = A * exp(l1+x[1][j]);
            u[i][0] = A * exp(x[0][i]);
            u[i][N2] = A * exp(l2+x[0][i]);
            u[0][j] = A * exp(x[1][j]);
            u[N1][j] = A * exp(l1+x[1][j]);
        }
    }
    int count = 0;
    double norm;
    do {
        count++;
        double** r = new double* [N1];
        for(int i(0); i<=N1; i++) {
            r[i] = new double[N2];
        }
        for(int i(1); i<N1; i++) {
            for(int j(1); j<N2; j++) {
                r[i][j] = Af(u, x, i, j) - phi;
            }
        }
        double scalR = 0;
        for(int i(1); i<N1; i++) {
            for(int j(1); j<N2; j++) {
                scalR += r[i][j] * r[i][j] * h1*h1;
            }
        }
        double scalA = 0;
        for(int i(1); i<N1; i++) {
            for(int j(1); j<N2; j++) {
                scalA += Af(r, x, i, j) * r[i][j] * h1*h1;
            }
        }
        double tau = scalR / scalA;

        //Y[i][j]
        for(int i(1); i<N1; i++) {
            for(int j(1); j<N2; j++) {
                y[i][j] = u[i][j] - r[i][j] * tau;
            }
        }

        //Norm
        norm = -1;
        for(int i(1); i<N1; i++) {
            for(int j(1); j<N2; j++) {
                if(fabs(y[i][j]-u[i][j]) > norm) {
                    norm = fabs(y[i][j]-u[i][j]);
                }
            }
        }

        //
        for(int i(1); i<N1; i++) {
            for(int j(1); j<N2; j++) {
                u[i][j]=y[i][j];
            }
        }
    } while(norm > Eps);

    for(int i(0); i<=N1; i++) {
        for(int j(0); j<=N2; j++) {
            cout<< setw(5) << setprecision(3) << func(x[0][i], x[1][j])<< " ";
        }
        cout<<endl<<endl;
    }
    cout << endl;
    cout<<"Count: " << count <<endl;
    for(int i(0); i<=N1; i++) {
        for(int j(0); j<=N2; j++) {
            cout<< setw(5) << setprecision(3) << y[i][j]<< " ";
        }
        cout<<endl<<endl;
    }
    system("pause");
    return 0;
}