#include <math.h>
#include <stdio.h>
#define MAXN 10000

/* data ulohy - globalni promenne */
double p = 0.25;
double a = 0;
double b = 1;
double T = 1;

int     n = 20;
int     maxk = 600;

double  h;
double  tau;

double  sigma;
 
double f(double x, double t);
double alfa(double t);
double beta(double t);
double u0(double x);

main()
{
    FILE *fid;
    double uk[MAXN], ukplus1[MAXN];
    int i, k;
    double xi, tk;

    h = (b - a) / n;
    tau = T / maxk;
    sigma = p * tau / (h * h);
    printf("sigma: %g\n", sigma);

    fid = fopen("reseni.dat", "w+");
    k = 0; 
    for(i = 0; i <= n; i++)
    {
        xi = a + i * h; uk[i] = u0(xi);
        fprintf(fid, "%18.8g ", uk[i]);
    }
    fprintf(fid, "\n");

    for(k = 0; k <= maxk; k++)
    {
        tk = k * tau;
        
        ukplus1[0] = alfa(tk);
        ukplus1[n] = beta(tk);

        for(i = 1; i < n; i++)
        {
            xi = a + i * h;
            ukplus1[i] = (1 - 2 * sigma) * uk[i] + sigma * uk[i - 1] + sigma * uk[i + 1] + tau * f(xi, tk);
        }
        for(i = 0; i <= n; i++)
        {
            uk[i] = ukplus1[i];
            fprintf(fid, "%18.8g ", ukplus1[i]);
        }
        fprintf(fid, "\n");
    }
    fclose(fid);
    printf("Reseni ulozeno do souboru: reseni.dat\n");
}

double f(double x, double t)
{
    return exp(-t);
}

double alfa(double t)
{
    return 0;
}

double beta(double t)
{
    return 0;
}

double u0(double x)
{
    return sin(M_PI * x);
}