Уравнение теплопроводности с переменным шагом

#include <iostream>
#include <fstream>
#include <string.h>  
#include <math.h>
#include <iomanip>
using namespace std;
int main()
{ 
    const int N = 20., lambda = 5.0;
    float gamma = 1.16126, dt = 0.01, dx = 0.001, E = 0.01;
    float u[N], u1[N], x[N], a[N], b[N], c[N], d[N];
    int i, r, t = 0;
 
    x[0] = 0.;
    for(i = 1; i <= N - 1; i++)
        x[i] = x[i - 1] + dx*pow(gamma,i - 1);
    FILE *F;
        char *file;
        file = "Results1.txt";
        F = fopen(file,"w");
        for(i = 0 ; i < N; i++)
            fprintf(F,"%f\n",x[i]);
 
    for(i = 0; i <= N - 2; i++)          //начальные и граничные условия:
        u[i] = u1[i] = 100.;      //при t = 0 и x=0      
    u[N - 1] = u1[N - 1] = 1000.;   //при x = N-1
    
    float z = 2*lambda*dt;
 
    while(true)
    {
        t = t + dt;
 
        for(i = 1; i <= N - 2; i++)
        {
            float hh = x[i + 1] - x[i - 1],
                h1 = x[i + 1] - x[i],
                h2 = x[i] - x[i - 1];
            u[i] = z/hh*(u1[i+1]/h1+u1[i-1]/h2)+u1[i]*(1-z/(h1*h2));
        }
        
        int j = 0;
        for(i = 1; i <= N - 2; i++)
        {
            r = abs((u[i]-u1[i])/u1[i]);
            if(r > j)
                j = r;
        }
        if(j > E)
        {
            u1[i] = u[i];
        }
        else
        {
            fprintf(F,"%f\n");
            for(i = 0 ; i < N; i++)
            {
                fprintf(F,"%f\n",u[i]);
                
            }
            break;
        }
    }
    
    fclose(F);
 
    return 0 ;
}

const int N = 20;
const double lambda = 5.0;

#include <iostream>
#include <fstream>
#include <cmath>
using namespace std;
int main()
{ 
    const int N = 20;
    double gamma = 1.16126;
    double x[N];
    x[0] = 0.;
    double pgamma=1.0,dx=0.001;
    for(int i = 1; i <= N - 1; i++)
    {
        x[i] = x[i - 1] + dx*pgamma;
        pgamma*=gamma;
    }
    char file[]="Results1.txt";
    ofstream F(file);
    if(!F)
        return 1;
    for(int i = 0 ; i < N; i++)
       F<<x[i]<<endl;
    double u[N],u1[N];
    for(int i = 0; i <= N - 2; i++)          //начальные и граничные условия:
        u[i] = u1[i] = 100.;      //при t = 0 и x=0      
    u[N - 1] = u1[N - 1] = 1000.;   //при x = N-1
 
    double lambda = 5.0;
    double dt = 0.01;
    double z = 2.0*lambda*dt;
    double t=0.;
    double eps=1e-4;
    while(true)
    {
        t += dt;
 
        for(int i = 1; i <= N - 2; i++)
        {
            double hh = x[i + 1] - x[i - 1];
            double h1 = x[i + 1] - x[i];
            double h2 = x[i] - x[i - 1];
            u[i] = z/hh*(u1[i+1]/h1+u1[i-1]/h2)+u1[i]*(1.-z/(h1*h2));
        }
        
        double j = 0;
        for(int i = 1; i <= N - 2; i++)
        {
            double r = fabs((u[i]-u1[i])/u1[i]);
            if(r > j)
                j = r;
        }
        if(j > eps)
        {
            for(int i = 0; i <= N - 1; i++)
            u1[i] = u[i];
        }
        else
        {
            F<<endl;
            for(int i = 0 ; i < N; i++)
               F<<u[i]<<endl;
            break;
        }
    }
    F.close();
    cout<<"Result in Results1.txt\n";
    system("pause");
    return 0 ;
}