Решение начально краевой задачи для уравнения теплопроводности

#include<iostream>
#include<istream>
#include<math.h>
#include<cmath>
#define PI 3.14159265358979323846
using namespace std;
const int N =10;
double f(double x)
{
    return sin(2 * PI*x) - 1 / (2 * PI*PI);
}
double phi(double x)
{
    return x*x;
}
double psi_0(double t)
{
    //t = 1;
    return -2 * PI*(1 - exp(-t));
}
double psi_1(double t)
{
    return 1;
}
 
int main()
{
    //double pi = 3.14159265358979323846;
    double h, tau, timing = 0, t = 0;
    double T = 1;
    h = 1.0 /(double) N;//посчитали шаг по x  
    tau =  h/ (double)4;//посчитали шаг по t  
    double Mas = T / tau; //сколько будет слоев, чтобы знать размерность массива  
    int l = 0;
    double B[N + 1];
    double u0[N + 1], u[N + 1];
    int z = 0; 
    double A[10*(N + 1)*(N + 1)][N + 1];
    
    for (int i = 0; i <= N; i++)
    {
        u0[i] = phi(i*h);
        A[0][i] = u0[i];
        cout << u0[i]<<endl;
    }
    //переопределяем значения на границах нулевого слоя  
    double a;// = (tau / (double)(4 * PI*h*h));
    u0[0] = psi_0(0);
    u0[N] = psi_1(0);
    cout << endl;
    while (timing <= T)
    {
        timing = timing + tau;
        z = z + 1;
        l = z;
        a=(timing/ (double)(4 * PI*h*h));
        for (int i = 1; i <= N-1; i++)
        {
            u[i] = (u0[i] + (a*(u0[i + 1] - 2 * u0[i] + u0[i - 1])))+ timing * f(i*h);
            A[l][i] = u[i];
        }
        u[0] = psi_0(timing);
        u[N] = -(u0[N] - u0[N - 1]) / (double)(h)+1;
        A[l][0] = u[0];
        A[l][N] = psi_1(N);//=u[N];//
 
        for (int i = 0; i <= N; i++)
        {
            u0[i] = u[i];
            cout << A[l][i] << " ";
        }
        cout << endl;
    }
 
    
    system("pause");
    return 0;
}