using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using System.IO;
namespace Laba11
{
    class Program
    {
        static double aa = 1.0, eps = 1e-10;
        static double xmin, xmax, tmin, tmax;
        static int n, m;
        static double h, t;
        static double[] x1;
        static double[] t1;
        static void Main(string
            [] args)
        {
            System.Threading.Thread.CurrentThread.CurrentCulture = new System.Globalization.CultureInfo("en-US");
            Console.Write("h: ");
            h = Convert.ToDouble(Console.ReadLine());
            Console.Write("t max: ");
            tmax = Convert.ToDouble(Console.ReadLine());
            Console.Write("t: ");
            t = Convert.ToDouble(Console.ReadLine());
            if (t - h * h / 2.0 < eps)
                Console.Write("Stable");
            else
                Console.Write("No stable");
            xmin = 0.0; xmax = 1.0; tmin = 0.0;
            _1();
            _2();
            Console.Read();
        }
        static void _1() //Вычисление явной схемы
        {
            StreamWriter wr = new StreamWriter("out1.txt");
            n = Convert.ToInt32((xmax - xmin) / h);
            m = Convert.ToInt32((tmax - tmin) / t);
            x1 = new double[n + 1];
            t1 = new double[m + 1];
            double[][] u = new double[n + 1][];
            for (int i = 0; i < n + 1; i++)
            {
                x1[i] = h * i;
                u[i] = new double[m + 1];
                u[i][0] = fi(x1[i]); // начальное условие
            }
            for (int j = 0; j < m + 1; j++)
                t1[j] = t * j;
            // вычисление значений искомой функции
            for (int j = 1; j < m + 1; j++)
            {
                u[0][j] = a1(t1[j]); // на левой границе
                for (int i = 1; i < n; i++)// центральное значение
                    u[i][j] = (
                        f(x1[i], t1[j - 1])
                        + aa * (u[i + 1][j - 1]
                            - 2.0 * u[i][j - 1]
                            + u[i - 1][j - 1]) / (h * h)
                        ) * t + u[i][j - 1];
                u[n][j] = a2(t1[j]); // на правой границе
            }
            // вывод значений
            for (int i = 0; i < n + 1; i++)
            {
                for (int j = 0; j < m + 1; j++)
                    wr.WriteLine(x1[i] + " " + t1[j] + " " + u[i][j]);
                wr.WriteLine();
            }
            wr.Close();
       }
        static void _2() //Вычисление неявной схемы
        {
            StreamWriter wr = new StreamWriter("out2.txt");
            n = Convert.ToInt32((xmax - xmin) / h);
            m = Convert.ToInt32((tmax - tmin) / t);
            x1 = new double[n + 1];
            t1 = new double[m + 1];
            double[] a = new double[n - 1];
            double[] b = new double[n - 1];
            double[] c = new double[n - 1];
            double[] d = new double[n - 1];
            double[] e = new double[n];
            double[] v = new double[n];
            double[][] u = new double[n + 1][];
            for (int i = 0; i < n + 1; i++)
            {
                x1[i] = h * i;
                u[i] = new double[m + 1];
                u[i][0] = fi(x1[i]); // начальное условие
            }
            for (int j = 0; j < m + 1; j++)
                t1[j] = t * j;
            // вычисление значений искомой функции (решение системы прогонкой)
            for (int j = 1; j < m + 1; j++)
            {
                u[0][j] = a1(t1[j]); // на левой границе
                u[n][j] = a2(t1[j]); // на правой границе
                a[0] = 0;
                b[0] = -h * h / t - 2.0 * aa;
                c[0] = aa;
                d[0] = -u[1][j - 1] * h * h / t - f(x1[1], t1[j]) * h * h - aa * u[0][j];
                for (int i = 1; i < n - 2; i++)// Нахождение коэфицентов
                {
                    a[i] = aa;
                    b[i] = -h * h / t - 2.0 * aa;
                    c[i] = aa;
                    d[i] = -u[i + 1][j - 1] * h * h / t - f(x1[i + 1], t1[j]) * h * h;
                }
                a[n - 2] = aa;
                b[n - 2] = -h * h / t - 2.0 * aa;
                c[n - 2] = 0;
                d[n - 2] = -u[n - 1][j - 1] * h * h / t - f(x1[n - 1], t1[j]) * h * h - aa * u[n][j];
                e[0] = v[0] = 0;
                for (int i = 0; i < n - 1; i++)// Нахождение промежуточных коэфицентов
                {
                    e[i + 1] = c[i] / (b[i] - a[i] * e[i]);
                    v[i + 1] = (a[i] * v[i] - d[i]) / (b[i] - a[i] * e[i]);
                }
                u[n - 1][j] = v[n - 1];
                for (int i = n - 1; i > 1; i--)// Вычисление значения искомой функции
                    u[i - 1][j] = e[i - 1] * u[i][j] + v[i - 1];
            }
            // вывод значений
            for (int i = 0; i < n + 1; i++)
            {
                for (int j = 0; j < m + 1; j++)
                    wr.WriteLine(x1[i] + " " + t1[j] + " " + u[i][j]);
                wr.WriteLine();
            }
            wr.Close();
        }
        static double f(double x, double t) // неоднородная часть
        {
            return x*x;
        }
        static double fi(double x) // начальное условие
        {
            return Math.Sin(Math.PI * x) * Math.Cos(x);
        }
        static double a1(double t) // левая граница
        {
            return 0;
        }
        static double a2(double t) // правая граница
        {
            return t;
        }
    }
}