程序代码：

//general scalar equation: d(rho*phi)/dt=d/dx(rho*D*(dphi/dx))
//W=West(i); C=centre(i+1); E=East(i+2)
//w=west(i+1/2); e=east(i+3/2)
#include<iostream>                                 
#include<fstream>                                      //for files
using namespace std;
void Square(int n, double dx, double a, double b, double x[], double y[])
{
     int i;
     for(i=0; i<n; i++)
    {
             x[i]=(i+0.5)*dx;                          //define initial x value for each cell, mid-point
             if (i>a&&i<b)                             //define initial y value for each cell 
             y[i]=1.0;                               
             else
             y[i]=0.0;                               
    }
}
double Diffusion(int i, double dx, double y[], double rho[], double D[])
{
     double yW, yC, yE, rhoW, rhoC, rhoE, dW, dC, dE, phiW, phiC, phiE, phiw, phie, rhodw, rhode;
     yW=y[i]; yC=y[i+1]; yE=y[i+2];                               //y at W, C and E
     rhoW=rho[i]; rhoC=rho[i+1]; rhoE=rho[i+2];                   //density at W, C and E
     dW=D[i]; dC=D[i+1]; dE=D[i+2];                               //diffusivity at W, C and E
     phiW=yW/rhoW; phiC=yC/rhoC; phiE=yE/rhoE;                    //y divided by density to get phi
     phiw=phiC-phiW; phie=phiE-phiC;                              //dleta phi
     rhodw=(rhoW*dW+rhoC*dC)/2; rhode=(rhoE*dE+rhoC*dC)/2;        //density*diffusivity at w and e
     return((rhode*phie-rhodw*phiw)/(dx*dx));                     //diffusion term
}
void Output(int n, double x[], double y[])
{
    int i;
    ofstream outdata;                                            //libray for output                   
    outdata.open("output.dat");                                  //open the file
    if(outdata.is_open())                                        //check if the file is open
    {
             for(i=0; i<n; i++)                                  //for each point
             {
                      outdata<<x[i]<<" "<<y[i]<<endl;            //write the value of x and y to file
             }
    outdata.close();                                             //close the file
    }
    else                                                         //if _open() return falsch
    {
             for(i=0; i<n; i++)
             {
                      cout<<x[i]<<" "<<y[i]<<endl;
             }
    }
}
int main()
{
    double l=1.0;                                      //total length
    double tpmax=1.0;                                  //total time
    double dt=0.01;                                    //time step interval
    const int n=100;                                   //define cell number
    double dx, a, b, tp;
    double x[n], y[n], rho[n], D[n], RHS[n], ybuffer[n];
    dx=l/n;                                            //calculate the length of each cell
    a=0.3*n;                                           //define the critical point
    b=0.7*n;                                           //define the critical point
    int i;
    for(i=0; i<n; i++)                                 //define the density and diffusivity values of air
    {
             rho[i]=1.225;
             D[i]=0.00002;
    }
                                                      
    Square(n, dx, a, b, &x[n], &y[n]);                               //switch on square
  
    for(tp=0.0; tp<tpmax; tp=tp+dt)
    {
              for(i=0; i<(n-1); i++)
              {
                       RHS[i+1]=Diffusion(i, dx, &y[n], &rho[n], &D[n]);                //switch on diffusion                               
                       ybuffer[i+1]=y[i+1]+dt*RHS[i+1];                                 //explicit time integration to get y
                       y[i+1]=ybuffer[i+1];
              }
    }
  
    Output(n, &x[n], &y[n]);                                         //results output
  
    return 0;                                                        // finish
}

Square(n, dx, a, b, x, y);