Составить таблицу значений функции y, заданной неявно - C (СИ)

Формулировка задачи:

Составить таблицу значений функции y, заданной неявно: x+y*y*y=e^(-x*x-y)+1, для 0<=x<=1 с шагом h=0,1 (написать на С, а не С++)

Решение задачи: «Составить таблицу значений функции y, заданной неявно»

Листинг программы

/*Составить таблицу значений функции y, заданной неявно:
x+y*y*y=e^(-x*x-y)+1, для 0<=x<=1 с шагом h=0,1
 
Преобразуем в
y = (exp(-y-x^2) + 1 - x)^(1/3)
х - параметр (константа на каждом отдельном шагу)
 
Решаем нелинейное ур-е относительно у методом итераций
Локализация корня: у Є [0;2]
dy = - exp(-y-x^2) / [3*(exp(-y-x^2) + 1 - x)^(2/3)]
При любом Х Є [0;1] и у Є [0;2] |dy| < 1
Более того, dy<0, т.е. процесс сходится с двух сторон, т.е. 
критерий сходимости: |Yk+1 - Yk| < e
*/
 
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
 
#define SEPARATOR "=================================================\n\n"
 
double func (double x, double y)
{   return pow((exp(-y-x*x) + 1 - x), 1/3.);
}
 
int main() 
{   double x = 0., dx = 0.1, e = 1.E-6;
    double ycur, ynext;
 
    printf("%10c%20c\n", 'X', 'Y');
    printf(SEPARATOR);
    while (x < 1.05)
    {   ycur = 0.;
        ynext = func (x, ycur);
        while (fabs (ynext - ycur) > e)
        {   ycur = ynext;
            ynext = func (x, ycur);
        }
        printf("%10.1lf%20.6lf\n", x, ynext);
        x += dx;
    }
    printf(SEPARATOR);
    getchar();
    return 0;
}

Объяснение кода листинга программы

В данном коде решается нелинейное уравнение методом итераций.
Неявная функция y задана формулой x+yyy=e^(-x*x-y)+1.
Для решения используется преобразование в y = (exp(-y-x^2) + 1 - x)^(1/3).
Для начала, итерации проходят в интервале [0;2] с шагом h=0.1.
Вычисляется функция y на каждом шаге и сравнивается с предыдущим значением y, пока разница не станет меньше заданной точности e=1.E-6.
Значения x и y выводятся на экран через каждые 0.1 единицы изменения x.
Значение x начинается с 0 и увеличивается до 1.05 с шагом 0.1.
Значение y на каждом шаге вычисляется как (exp(-y-x^2) + 1 - x)^(1/3).
Корень локализован в интервале [0;2].
При любом x в интервале [0;1] и любом y в интервале [0;2], |dy| < 1, что гарантирует сходимость процесса.
Критерий сходимости: |Yk+1 - Yk| < e.