restart;

Odegenerator := proc(V, x, y, df, const_values)
    local input_args, xi, F, result, a0, a1, a2, b0, b1, sol, Fsol, rows, numrows, eq, count, odeplot_cmd, ode_type, row_number, values;
    uses plots, PDEtools;
       if nargs = 1 and V = "help" then
        printf("Use this procedure as follows:\n");
        printf("Define an ODE template:\n");
        printf("Odegenerator(V, x, y, df, const_values)\n");
        printf("V: A set of values for iteration over constants (if df > 0)\n");
        printf("x: The independent variable\n");
        printf("y: The function\n");
        printf("df: The row number in the DataFrame or 0 for manual input\n");
        printf("const_values: A list of values for the constants (used if df = 0)\n");
        return;
    end if;

    if nargs < 4 or nargs > 5 then
        error "Incorrect number of arguments. Expected: V, x, y, df, [const_values (optional)]";
    end if;

    # Determine the ODE type using odeadvisor for the global eq_template
    ode_type := odeadvisor(eq_template);

    # Display the ODE and its type
    print(eq_template, ode_type);

    rows := [];
    count := 0;
    ###################### BOF manuele invoer ###################
    if df = 0 then
    # If df = 0, use const_values for substitution
    if nargs < 5 or not type(const_values, list) then
        error "When df = 0, a list of constant values must be provided as the fifth argument.";
    end if;

    # Assign constant values
    if nops(const_values) <> 5 then
        error "The list of constant values must contain exactly 5 elements.";
    end if;

    # Find the corresponding row number by unique identification
    count := 1;
    for a0 in V do
        for a1 in V do
            for a2 in V do
                for b0 in V do
                    for b1 in V do
                        if [a0, a1, a2, b0, b1] = const_values then
                            row_number := sprintf("%d", count);  # Convert to string
                        end if;
                        count := count + 1;
                    end do;
                end do;
            end do;
        end do;
    end do;

    if not assigned(row_number) then
        row_number := "Unique (outside iterative rows)";  # Mark as unique
    end if;

    # Substitute the given values
    eq := subs({'a__0' = const_values[1], 'a__1' = const_values[2], 'a__2' = const_values[3], 'b__0' = const_values[4], 'b__1' = const_values[5]}, eq_template);

    # Solve the equation
    sol := dsolve(eq, y(x));
    if type(sol, `=`) then
        Fsol := rhs(sol);
    else
        Fsol := "No explicit solution";
    end if;

    # Display the solution and its row number
    odeplot_cmd := DEtools[DEplot](eq, y(x), x = 0 .. 2, y = -10 .. 10, [[y(0) = 1]]);
    print(plots:-display(odeplot_cmd, size = [550, 550]));

    printf("The found function is:\n");
    print(Fsol);
    printf("The corresponding row number is: %s\n", row_number);

    # -- Start of Additional Functionality --
    # Optionally display the simplified ODE
    printf("The simplified ODE using the given coefficients is:\n");
    print(eq, ode_type);
    # -- End of Additional Functionality --

    return Fsol;
     ################# EOF manuele berekening ##################
     ############## BOF iterative berekening##############
    else
        # Iterative approach for DataFrame generation
        for a0 in V do
            for a1 in V do
                for a2 in V do
                    for b0 in V do
                        for b1 in V do
                            xi := x;
                            F := y;

                            # Substitute constant values into the ODE
                            eq := subs({'a__0' = a0, 'a__1' = a1, 'a__2' = a2, 'b__0' = b0, 'b__1' = b1}, eq_template);

                            sol := dsolve(eq, F(xi));
                            if type(sol, `=`) then
                                Fsol := rhs(sol);
                            else
                                Fsol := "No explicit solution";
                            end if;

                            rows := [op(rows), [a0, a1, a2, b0, b1, Fsol]];
                        end do;
                    end do;
                end do;
            end do;
        end do;

        numrows := nops(rows);
        result := DataFrame(Matrix(numrows, 6, rows), columns = ['a__0', 'a__1', 'a__2', 'b__0', 'b__1', y(x)]);

        interface(rtablesize = numrows + 10);

        if df > 0 and df <= numrows then
            a0 := result[df, 'a__0'];
            a1 := result[df, 'a__1'];
            a2 := result[df, 'a__2'];
            b0 := result[df, 'b__0'];
            b1 := result[df, 'b__1'];

            eq := subs({'a__0' = a0, 'a__1' = a1, 'a__2' = a2, 'b__0' = b0, 'b__1' = b1}, eq_template);

            # Display the additional parameters
            print(eq, ode_type, [df], [a0, a1, a2, b0, b1]);

            # Retrieve the solution
            Fsol := result[df, y(x)];

            # Display the solution in DEplot
            odeplot_cmd := DEtools[DEplot](eq, y(x), x = 0 .. 2, y = -10 .. 10, [[y(0) = 1]]);
            print(plots:-display(odeplot_cmd, size = [550, 550]));

            printf("The found function for row number %d is:\n", df);
            print(Fsol);

        else
            printf("The specified row (%d) is out of bounds for the DataFrame.\n", df);
        end if;

        return result;
     ########## EOF iteratief bwrekening ########################
    end if;

end proc:

# Test cases
V := {0, 1};
eq_template := diff(y(t), t) = 'a__0'*sin(t) + 'a__1'*y(t) + 'a__2'*y(t)^2 + 'b__0'*exp(-t);