1次元拡散方程式の陽解法(C#のコードと結果)

以下工事中です。

using System;
using System.IO;
using System.Text;

class Program
{
    static void Main()
    {
        int    n_max, j_max;
        double dt, dx, t_max;
        double r = 0.25;     //diffusion number

        //シミュレーションの説明
        Console.WriteLine("This is a simulation of the one-dimensional diffusion equation.");
        Console.WriteLine("The Dirichlet boundary condition is imposed at the spatial boundary x=0 and x=1.");
        Console.WriteLine(string.Format("To stabilize the solution, we fix the diffusion number to {0}.", r));

        //シミュレーションパラメータのコンソール入力
        Console.Write("Input a number of time step: ");
        n_max = int.Parse(Console.ReadLine());
        Console.Write("Input a number of spatial step: ");
        j_max = int.Parse(Console.ReadLine());

        //空間と時間の解像度
        dx = 1.0 / j_max;
        dt = r * Math.Pow(dx,2); //diffusion number の定義からdtが定まる
        t_max = n_max * dt;
        Console.WriteLine(string.Format("The simulation time is {0}", t_max));

        //格子点の作成
        double[] t = new double[n_max + 1];
        double[] x = new double[j_max + 1];
        for(int n = 0; n <= n_max; n++)
        {
            t[n] = n * dt;
        }
        for(int j = 0; j <= j_max; j++)
        {
            x[j] = j * dx;
        }
        
        //格子点上の場の生成(uの引数は空間格子番号、uをtime step毎に更新していく)
        double[] u = new double[j_max + 1];
        //初期条件の代入
        for(int j=0; j<=j_max; j++)
        {
            u[j] = 2.0 * x[j] * (1.0 - x[j]);
        }

        //csvファイルインスタンス生成
        StreamWriter file = new StreamWriter(@"/Users/kaishu/Projects/diffusion_explicit_ver2/data.csv", false, Encoding.UTF8);
        //カラム名の書き込み
        file.Write("Label");
        file.Write(",");
        file.Write("Time");
        file.Write(",");
        for (int j = 0; j <= j_max; j++)
        {
            if(j<=j_max - 1)
            {
                file.Write(string.Format("u_" + "{0}", j));
                file.Write(",");
            }
            else
            {
                file.WriteLine(string.Format("u_" + "{0}", j));
            }
        }

        //陽解法による数値積分
        for(int n =0; n<=n_max; n++)
        {
            //ラベルと時間の記録
            file.Write(string.Format("{0}", n));
            file.Write(",");
            file.Write(string.Format("{0}", t[n]));
            file.Write(",");

            //境界条件
            u[0] = 0;
            u[j_max] = 0;

            file.Write(string.Format("{0}", u[0]));
            file.Write(",");

            if(n==0)
            {
                for(int j = 1; j <= j_max - 1; j++) 
                {
                    file.Write(string.Format("{0}", u[j]));
                    file.Write(",");
                }
            }
            else
            {
                for(int j = 1; j <= j_max-1; j++)
                {
                    u[j] = r * u[j + 1] + (1.0 - 2.0 * r) * u[j] + r * u[j - 1];
                    file.Write(string.Format("{0}", u[j]));
                    file.Write(",");
                }
            }
            file.WriteLine(string.Format("{0}", u[j_max]));
        }

        file.Close();
    } 
}