ここでは初期条件として温度が100、境界条件としてまわりの温度が0である2次元の系を考える。

#include <stdio.h>

#define N 100           //正方形である2次元の系のそれぞれの辺の長さに値する

int main()
{
        int t,j,k,tmax;
        double T[2][N][N];

        FILE *output;
        output=fopen("output.data","w");

/*最大繰り返し回数とは、時間的にどれくらい経過させるかのこと。*/
        printf("最大繰り返し回数を入力してください。\n");
        scanf("%d",&tmax);

/*初期条件として全体の温度100*/
        for(j=0;j<N;j++) {
                for(k=0;k<N;k++) {
                        T[0][j][k]=100.0;
                }
        }

        for(t=0;t<tmax;t++) {

/*境界条件として回りの温度0*/
                for(j=0;j<N;j++) {
                        T[0][0][j]=0.0;
                        T[0][N-1][j]=0.0;
                        T[0][j][0]=0.0;
                        T[0][j][N-1]=0.0;
                }

                for(j=1;j<N-1;j++) {
                        for(k=1;k<N-1;k++) {
                                T[1][j][k]=T[0][j][k]+0.25*(T[0][j+1][k]+T[0][j-1][k]+T[0][j][k+1]+T[0][j][k-1]-4.0*T[0][j][k]);
                        }
                }

/*時間的に過去の情報は破棄するためと、この時間的forループの計算のために配列Tのデータ処理*/
                for(j=1;j<N-1;j++) {
                        for(k=0;k<N-1;k++) {
                                T[0][j][k]=T[1][j][k];
                        }
                }

        }

/*最終結果には境界条件である温度0の情報がはいってないため、ここで計算する。*/
        for(j=0;j<N;j++) {
                T[1][0][j]=0.0;
                T[1][N-1][j]=0.0;
                T[1][j][0]=0.0;
                T[1][j][N-1]=0.0;
        }

        for(j=0;j<N;j++) {
                for(k=0;k<N;k++) {
                        fprintf(output,"%d %d %f\n",j,k,T[1][j][k]);
                }
                fprintf(output,"\n");
        }

        fclose(output);

        return 0;
}

上のプログラムを実行して、最大繰り返し回数を1000と入力した場合のGNUPLOT出力結果である。

上の図はそれぞれの横軸がxy軸であり、縦軸が温度の高さを表す。下の図はxy平面を温度の高さによって、色を塗り分けたものである。