对这个抛物方程用古典差分格式，下面是源代码。但是有个问题让我很困惑。。。就是下面输出时要是没有if语句，结果很奇怪。
#include<stdio.h>
#include<math.h>
#define M 10    /*将空间区域等分为M份，即空间步长为1.0/M,有M+1个节点*/
#define N 100 /*将时间区域等分为N份，共N+1个时间节点；*/
void main()
{
    int i,j,k,m=0;
    double a[M-1][M-1]={0},x[M-1]={0},b[M-1]={0},h=1.0/M,t=1.0/N,r;
    r=t/h/h;
    for(i=1;i<M-2;i++){
        a[i][i-1]=r;a[i][i]=1-2*r;a[i][i+1]=r;
                          }
        a[0][0]=1-2*r;a[0][1]=r;
        a[M-2][M-3]=r;a[M-2][M-2]=1-2*r;
        for(i=0;i<M-1;i++)x[i]=exp(1.0*(i+1)/M);
        
         for(k=1;k<N+1;k++){  /*用这个循环来控制迭代，由于时间步长为1/200,即等分为200段，有201个节点，因此一共进行了201次迭代；*/
            double c[M-1]={0};/*下面构造关于边值的向量*/
            c[0]=r*exp((k-1)*t);
                        c[M-2]=r*exp(1+(k-1)*t);
  for(i=0;i<M-1;i++){  /*这个循环是用来进行对第k时间层的数值解向量的输出；因是先进行迭代运算再输出，因此本程序中没有输出初值向量，即第0层的数值解*/
                for(j=0;j<M-1;j++)c[i]+=a[i][j]*x[j];/*此循环用来求出第k时间层的第i个分量*/
            
                b[i]=c[i];/*用b[i]保存得到的c[i],由于在进行最外层循环时需要x[i],因此将b[i]保存新的x[i]*/
                m++;/*根据m最后的值可以判断一共进行了多少次循环；*/
if(k>=1&&k<=18&&i==4)printf("%22.17lf\n",b[i]);/*要是没有if控制语句，输出的结果非常奇怪。*为什么结果会有很多的不同*？/
            
                }
        for(i=0;i<M-1;i++)x[i]=b[i];/*这个赋值运算得到了新的x[i]，于是便可进行下一次的迭代*/
                
            }
        printf("%d\n",m);
        }