Derivatives: tables and graphs

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define DO 1
#define NOT 0
// диапазон интегрирования
#define FROM 0
#define TO 1
// настройки графика, отражающие масштаб по осям
#define Xscale 0.02
#define Yscale 0.02

float f(float x) { return x*x*cos(x); }

int cases; // количество строк, при вычислении с шагом h
float M0, M1, M2; // максимальные по модулю величины для производных соответствующих порядков

/* график табличной функции table[], гор. - table[0], верт table[colomn]*/
void draw(float *xarray, float *array, float Xstart, float Xend, float Ystart, float Yend){
    char **graph;
    int i, j;

    const int Xsize=(Xend-Xstart)/Xscale+1; //вычисление размера графика
    const int Ysize=(Yend-Ystart)/Yscale+1; // плюс два символа: на концец строки и рамки

    /* Выделяем память, заполняем поле пробелами и добавляем рамки */
    graph = malloc(sizeof(char*)*Ysize);
    for(i=0; i<Ysize; i++){
        char filler;
        graph[i] = malloc(sizeof(char)*Xsize);
            graph[i][Xsize] = 0; //помещаем в конец строки нулевой символ
            if(i) // если строка не последняя
                filler = ' '; // то просто заполняем поле
            else
                filler = '-'; //то рисуем горизонтальную линию
            graph[i][0]='|';
            for(j=1; j<Xsize; j++)
                graph[i][j]=filler; //заполняем пробелами пустые места;
    }
    /* отмечаем значения функций */
    // graph[строка - значение функции][столбец - аргумент]
    for(i=0; i<cases; i++)
        graph[ (int)(array[i]/Yscale) ][ (int)(xarray[i]/Xscale) ]='*'; // точка значения
    /* выводим график и сразу же освобождаем память */
    for(i=Ysize-1; i>=0; i--){
        printf("%.3f %s\n", i*Yscale, graph[i]);
        free(graph[i]);
    }
    free(graph);
    // вывести масштаб по горизонтали
    printf("      %d", FROM);
    for(i=1; i<Xsize-1; i++)
        printf(" ");
    printf("%d\n", TO);
}

/* h - шаг вычисления значений, int - режим работы (0 - без вывода таблицы) */
float derivative(float h, int mode){
    float *x, *y, *yL, *yR, *yC, *y2;

    if(mode)
        printf("x \t f(x) \t f'L(x) \t f'R(x) \t f'C(x) \t f''C(x)\n");
    cases = (TO-FROM)/h;
    int i;
    float x_;
    // значения
        x = malloc(sizeof(float)*cases);
        y = malloc(sizeof(float)*cases);
        yL = malloc(sizeof(float)*cases);
        yR = malloc(sizeof(float)*cases);
        yC = malloc(sizeof(float)*cases);
        y2 = malloc(sizeof(float)*cases);
    for(i=0, x_=FROM; i<=cases; x_+=h, i++){
        y[i] = f( x[i] = x_ ); // значение аргумента и функции
        if(y[i] > M0)
            M0 = y[i]; // выбираем максимальное значение функции
    }
    /* значения первой и второй производных */
    if(mode)
            printf("%.3f \t %.3f \t NaN    \t %.4f \t NaN    \t NaN\n",
                    x[0],   y[0],             yR[0]=(y[1]-y[0])/h);
    for(i=1; i<cases; i++){
        /* знач. первой производной (левая аппрокс.) */
        yL[i] = (y[i]-y[i-1])/h;
        /* знач. первой производной (правая аппрокс.) */
        yR[i] = (y[i+1]-y[i])/h;
        /* знач. первой производной (центральная аппрокс.) */
        yC[i] = (y[i+1]-y[i-1])/(2*h);
        if(yC[i] > M1)
            M1 = yC[i]; // выбираем максимальное значение M1
        /* знач. второй производной */
        y2[i] = (y[i-1]-2*y[i]+y[i+1])/(h*h);
        if(y2[i] > M2)
            M2 = y2[i]; // выбираем максимальное значение M2
        /* вывод таблицы */
        if(mode)
            printf("%.3f \t %.3f \t %.4f \t %.4f \t %.4f \t %.3f\n",
                x[i], y[i], yL[i], yR[i], yC[i], y2[i]);
    }
    if(mode){
        /* График функции и график первой производной */
        printf("График функции: ордината - f(x), абсцисса - x\n");
        draw(x, y, FROM, TO, 0, M0);
        printf("График первой производной функции: ордината - f'(x), абсцисса - x\n");
        draw(x, yL, FROM, TO, 0, M1);
    }
    return yC[1];
}

short int Runge(double eps, float h){
    // считаем производные
    float val1 = derivative(h, NOT);
    float val2 = derivative(h/2, NOT);
    // если разница мала то 1, иначе 0
    if( fabsf(val1 - val2) < eps ){
        printf("Шаг %f подошёл\n", h);
        return 1;
    }else{
        return 0;
    }
}

int main(){
    float h, eps;
    srand(time(NULL)); //let's be random

    printf("Введите начальный шаг: ");
    scanf("%f", &h);
    printf("Точность для принципа Рунге: ");
    scanf("%f", &eps);

    /* определение максимальных значений */
    printf("Определение максимальных значений...\n");
    float temp = derivative(h, 0);
    printf("Максимальные значения определены: M1=%f, M2=%f\n", M1, M2);
    /* запуск метода Рунге */
        printf("Запуск метода Рунге.\n");
    while( Runge(eps, h) == 0 ) h/=2;
    // рассчитываем производные в этой точке и строим графики
    derivative(h, DO);

    return 0;
}