chenged proga

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<iostream>
#include<vector>
#include<string>
#include<sstream>
#include<fstream>

using namespace std;

 //const double alpha = 0;
 const double N = 20;
 const double h = 1.0 / N;
 const double PI = 3.141592653589793238463;
 const double eps = 0.01;
 //const double tau = h * h / sin(PI * h);
 const double tau = h;

 double alpha(double x, double y) {
     //return 0;
     //return x;
     if(x == 0) return sin(PI * y);
     else if(y == 0) return sin(PI * x);
     else return 0;

 }
 double f(double x, double y)
 {
     return x * y;
     //return 2 * PI * PI * sin(PI * x) * sin(PI * y);
     //return x*y;

     //return 2 * PI * sin(PI * x) * sin(PI * y) - y * PI * cos(PI * x) * sin(PI * y) - x * PI * cos(PI * y) * sin(PI * x);
 }
 double f(int i, int j)
 {
     return f((double)i * h, (double)j * h);
 }
 double Real_U(double x, double y)
 {
     return sin(PI * x) * sin(PI * y);
     //return sin(PI * x) * sin(PI * y) * x * y;
 }
 double Real_U(int i, int j)
 {
     return Real_U((double)i * h, (double)j * h);
 }

 int indic(int i, int j) {
     if(i == 0 || i == N || j == 0 || j == N) return 1;
     else return 0;
 }
 string to_string(double a){
     ostringstream strs;
     strs << a;
     string str = strs.str();
     return str;
 }
 void PrintToFile(string filename,string str){
    const   char* s = filename.c_str();
     std::ofstream out(s);
    out << str;
    out.close();
 }

 vector<double> ThomasAlgoritm(vector<double> A, vector<double> C, vector<double> B, vector<double> F)
 {
     vector<double> C_;   C_.push_back(C[0]);
     vector<double> F_;   F_.push_back(F[0]);
     vector<double> item;   item.push_back(0.0);
     int n = C.size();
     for (int i = 1; i < n; i++)
     {
         C_.push_back(C[i] - A[i - 1] * B[i - 1] / C_[i - 1]);
         F_.push_back(F[i] - A[i - 1] * F_[i - 1] / C_[i - 1]);
         item.push_back(0.0);
     }
     item[n - 1] = F_[n - 1] / C_[n - 1];
     for (int i = n - 2; i >= 0; i--)
     {
         item[i] = (F_[i] - B[i] * item[i + 1]) / C_[i];
     }
     return item;
 }
 class approx
 {
     vector<vector<double> > sh_k;
     vector<vector<double> > sh_k_half;

 public:
     approx()
     {
         sh_k = *(new vector<vector<double> >);
         sh_k_half = *(new vector<vector<double> >);

         int i, j;
         for (i = 0; i <= N; i++)
         {
             sh_k.push_back(*(new vector<double>));
             sh_k_half.push_back(*(new vector<double>));
             for (j = 0; j <= N; j++)
             {
                 sh_k[i].push_back(alpha(i * h, j * h));
                 sh_k_half[i].push_back(alpha(i * h, j * h));
             }
         }
     }
     void step_2(){
         int i, j;
         vector<double> solve = *(new vector<double>);
         vector<double> A = *(new vector<double>);
         vector<double> C = *(new vector<double>); C.push_back(1.0);
         vector<double> B = *(new vector<double>); B.push_back(0.0);
         vector<double> F = *(new vector<double>);
         for (i = 0; i < N - 1; i++)
         {
             A.push_back(tau / (2 * h * h));
             B.push_back(tau / (2 * h * h));
             C.push_back(-(tau / (h * h) + 1));
         }
         A.push_back(0.0); C.push_back(1.0);
         for (j = 1; j < N; j++)
         {
             // заполняем F
             F.clear();
             F.push_back(alpha(0, j * h));
             //for (i = 1; i < N; i++) F.push_back(alpha(i * h, j * h));
             for (i = 1; i < N; i++)
             {
                 F.push_back(-sh_k[i][j] - tau * ((sh_k[i][j + 1] - 2 * sh_k[i][j] + sh_k[i][j - 1]) / (h * h) + f(i, j)) / 2);
             }
             F.push_back(alpha(1, j * h));
             // решаем систему
             solve = ThomasAlgoritm(A, C, B, F);
             for (i = 0; i <= N; i++)
                 sh_k_half[i][j] = solve[i];
         }
     }
     void step_3()
     {
         int i, j;
         vector<double> solve = *(new vector<double>);
         vector<double> A = *(new vector<double>);
         vector<double> C = *(new vector<double>); C.push_back(1.0);
         vector<double> B = *(new vector<double>); B.push_back(0.0);
         vector<double> F = *(new vector<double>);
         for (i = 0; i < N - 1; i++)
         {
             A.push_back(tau / (2 * h * h));
             B.push_back(tau / (2 * h * h));
             C.push_back(-(tau / (h * h) + 1));
         }
         A.push_back(0.0); C.push_back(1.0);
         for (i = 1; i < N; i++)
         {
             // заполняем F
             F.clear();
             F.push_back(alpha(i * h, 0));
             for (j = 1; j < N; j++)
             {
                 F.push_back(-sh_k_half[i][j] - tau * ((sh_k_half[i + 1][j] - 2 * sh_k_half[i][j] + sh_k_half[i - 1][j]) / (h * h) + f(i, j)) / 2);
             }
             F.push_back(alpha(i * h, 1));
             // решаем систему
             solve = ThomasAlgoritm(A, C, B, F);
             for (j = 0; j <= N; j++)
                 sh_k[i][j] = solve[j];
         }
     }
     void step_4()
     {
         step_2();
         step_3();
     }


     string GiveWolframString()
     {
         string s = "";
         int i, j;
         s += "ListPlot3D[{";
         for (i = 0; i <= N; i++)
             for (j = 0; j <= N; j++)
             {
                 s += "{";
                 s += to_string(i * h);
                 s += ",";
                 s += to_string(j * h);
                 s += ",";
                 s += to_string(sh_k[i][j]);
                 if(i + j < 2 * N) s += "}, ";
                 else s += "}}]";
             }
         return s;
     }
     string GiveWolframErrorString()
     {
         string s = "";
         int i, j;
         for (i = 0; i <= N; i++)
             for (j = 0; j <= N; j++)
             {
                 s += "{";
                 s += to_string(i * h);
                 s += ",";
                 s += to_string(j * h);
                 s += ",";
                 s += to_string(abs(Real_U(i, j) - sh_k[i][j]));
                 s += "}, ";
             }
         return s;
     }
     double error(){
         int i, j;
         double er = 0.0;
         for (i = 0; i <= N; i++)
             for (j = 0; j <= N; j++)
                 er += abs(sh_k[i][j] - Real_U(i, j)) * abs(sh_k[i][j] - Real_U(i, j));
         return sqrt(er)/N;
     }

 };

int main(void)
{
    class approx chema;
    double er = chema.error();
    int i;



    for ( i = 0; i < 200; i++)
    {
        chema.step_4();
        if (abs(er - chema.error()) < eps) break;
        //else
        er = chema.error();
        cout << i << "["<< er<< "] , ";
        //Console.WriteLine("{0} - {1}",i, er);
    }

    chema.step_4();

    string wolf_str = chema.GiveWolframString();
    string wolf_err = chema.GiveWolframErrorString();
    string wolf_item = "{" + to_string(log(N)) + "," + to_string(-log(er)) + "}";

    PrintToFile("wolf_str.txt",wolf_str);
    PrintToFile("wolf_err.txt",wolf_err);
    PrintToFile("wolf_item.txt",wolf_item);
}