Численное решение задачи Коши, результат выводится неправильно - C (СИ)

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

не понимаю в чем проблема задание было такое Построить алгоритмы численного решения задачи Коши явным, неявным методом Эйлера и методом Рунге-Кутта с автоматическим выбором шага интегрирования для достижения заданной локальной относительной погрешности  либо локальной абсолютной погрешности . Апостериорную оценку погрешности для явного метода Эйлера провести по правилу Рунге либо по оценке производной, входящей в выражение для погрешности. 2. Составить рабочую программу с использованием универсальных функций численного решения задачи Коши. Минимальный набор параметров функции решения задачи Коши должен включать: имя функции, вычисляющей производную неизвестной функции (правую часть дифференциального уравнения); интервал времени, на котором нужно найти решение; начальное условие; погрешность решения; массив временных отсчетов, в которых найдено решение; массив решений; максимальное число временных отсчетов; количество фактически полученных временных отсчетов, в которых найдены решения с заданной погрешностью.

                    Листинг программы
                
 
                            
                        
#define _USE_MATH_DEFINES
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
int a1;
 
float TEC(float t, float U)
{
    float z, e, E = 5, id, i0 = 3e-9, f = 1000, R = 5000, c = 1e-6, m = 1.7, fi = 0.026;
    a1++;
    e = E*sin(2 * (M_PI)*f*t + (M_PI) / 4);
    id = i0*(exp((e - U) / (fi*m)) - 1);
    z = id / c - U / (R*c);
    return z;
 
}
int RungeKutt(float f(float x, float y), float *t, float *y, float A, float B, float ny, float E, int Nmax)
{
    int i, j, n = 0;
    float q = 2, R, y1, y2, y3, t1, t2, t3, h, a[6] = { 0, (1. / 4), (3. / 8), (12. / 13), 1., (1. / 2) }, c[6] = { 16. / 135, 0, 6656. / 12825, 28561. / 56430, -9. / 50, 2. / 55 }, c1[6] = { 25. / 216, 0, 1408. / 2565, 2197. / 4104, -1. / 5, 0 }, b[5][5] = { { 1. / 4, 0, 0, 0, 0 }, { 3. / 32, 9. / 32, 0, 0, 0 }, { 1932. / 2197, -7200. / 2197, 7296. / 2197, 0, 0 }, { 439. / 216, -8., 3680. / 513, -845. / 4104, 0 }, { -8. / 27, 2., -3544. / 2565., 1859. / 4104, -11. / 40 } }, k[6];
    h = B - A;
    t1 = A;
    y1 = ny;
    do
    {
        y3 = y1;
        t3 = t1;
        for (i = 0; i<6; i++)
        {
            t2 = t1;
            y2 = y1;
            t2 += a[i] * h;
            for (j = 0; j<i; j++)
            {
                y2 += h*b[i - 1][j] * k[j];
            }
            k[i] = f(t2, y2);
        }
        for (i = 0; i<6; i++)
            y1 += h*c[i] * k[i];
        t1 += h;
        R = 0;
        for (i = 0; i<6; i++)
            R += (c[i] - c1[i])*k[i];
        R = h*R;
        if (fabs(R) <= E*fabs(y1))
        {
            y[n] = y1;
            t[n] = t1;
            if (fabs(R)<(E*fabs(y1) / 64))
                h = h*q;
            if ((t1 + h)>B)
                h = B - t1;
            n++;
        }
        else
        {
            h = h / q;
            y1 = y3;
            t1 = t3;
        }
    } while ((t1<B) && (n<Nmax));
    return n;
}
int yavniyEuler(float f(float x, float z), float *t, float *y, float a, float b, float ny, float E, float D, int Nmax)
{
    int i, n = 0, q = 2;
    float h, y1, t1, y2, y3, ny1, nt1, d, k = 0, R;
    h = b - a;
    t1 = a;
    y1 = ny;
    do
    {
        ny1 = y1;
        nt1 = t1;
        if (k == 0)
        {
            d = f(t1, y1);
            y2 = y1 + h*d;
        }
        k = 1;
        for (i = 0; i<q; i++)
        {
            if (i == 0)
            {
                y1 = y1 + (h / q)*d;
                y3 = y1;
            }
            else
                y1 = y1 + (h / q)*f(t1, y1);
            t1 += (h / q);
        }
        R = E*fabs(y1);
        if (y2 == 0)
        {
            R = R + D;
        }
        if (fabs(y1 - y2) <= (R))
        {
            k = 0;
            y[n] = y1;
            t[n] = t1;
            if (fabs(y1 - y2)<E*fabs(y1) / 4)
                h = q*h;
            n++;
        }
        else
        {
            h = h / q;
            y2 = y3;
            y1 = ny1;
            t1 = nt1;
        }
    } while (t1<b && n<Nmax);
    return (n);
}
int neyavniyEuler(float f(float x, float z), float *t, float *y, float a, float b, float ny, float E, int Nmax)
{
    int n = 0, q = 2;
    float h, h1, d, em = 0.000000119, y1, t1, y2, t2, y3, y4, f1, f2, v, R, R1;
    h = b - a;
    y1 = ny;
    t1 = a;
    y2 = a;
    d = f(a, ny);
    do
    {
        y4 = y1 + h*d;
        t2 = t1 + h;
        y3 = y2;
        if (y2 != 0)
            h1 = y2*sqrt(em);
        else
            h1 = sqrt(em);
        f1 = y2 - y1 - h*f(t2, y2);
        f2 = (y2 + h1) - y1 - h*f(t2, (y2 + h1));
        y2 = y2 - f1*h1 / (f2 - f1);
        while (fabs(y2 - y3)>E*fabs(y2))
        {
            f2 = f1;
            f1 = y2 - y1 - h*f(t2, y2);
            v = y2;
            y2 = y2 - (y2 - y3)*f1 / (f1 - f2);
            y3 = v;
        }
        R = (y4 - y2) / 2;
        if (fabs(R) <= E*fabs(y2))
        {
            y1 = y2;
            t1 = t2;
            y[n] = y2;
            t[n] = t2;
            d = f(t1, y1);
            if (fabs(R)<E*fabs(y2) / 4)
                h = h*q;
            if ((t2 + h)>b)
                h = b - t2;
            n++;
        }
        else
        {
            h = h / q;
        }
    } while (t1<b && n<Nmax);
    return n;
}
int main(void)
{
    int j, i, Nmax = 10000;
    float *t, *y;
    t = malloc(Nmax*sizeof(float));
    y = malloc(Nmax*sizeof(float));
    a1 = 0;
    i = yavniyEuler(TEC, t, y, 0, 1.5e-3, 0, 1e-6, 0.01, 10000);
    printf("%.8f n=%d a1=%d\n\n", y[i - 1], i, a1);
    a1 = 0;
    i = neyavniyEuler(TEC, t, y, 0, 1.5e-3, 0, 1e-6, 0.01, 10000);
    printf("%.8f n=%d a1=%d\n\n", y[i - 1], i, a1);
    a1 = 0;
    i = RungeKutt(TEC, t, y, 0, 1.5e-3, 0, 1e-6, 10000);
    printf("%.8f n=%d a1=%d\n\n", y[i - 1], i, a1);
}

Решение задачи: «Численное решение задачи Коши, результат выводится неправильно»

textual
                            Листинг программы 
                        
float f(float x, float y)