From 7853dcb74b414d1b0199a65592a289e09a0e916a Mon Sep 17 00:00:00 2001
From: Gabriel <sondar@50ND4R-NOTE.localdomain>
Date: Wed, 17 May 2023 11:13:02 +0300
Subject: [PATCH 01/24] comentario

---
 Energy_Preserving.jl | 500 +++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 500 insertions(+)
 create mode 100644 Energy_Preserving.jl

diff --git a/Energy_Preserving.jl b/Energy_Preserving.jl
new file mode 100644
index 00000000..4bc481c0
--- /dev/null
+++ b/Energy_Preserving.jl
@@ -0,0 +1,500 @@
+# Load the packages we will use.
+# These must first be installed using: import Pkg; Pkg.add("package_name")
+#import Pkg; Pkg.add("BSeries")
+using BSeries
+using Latexify  # Only needed for some pretty-printing cells below using `latexify`
+import SymPy; sp=SymPy;
+using LightGraphs
+using Combinatorics
+using RootedTrees
+using LinearAlgebra
+
+#Parametros iniciales
+#stisfies for order <5, and not for 5.
+A = [ 0 0 0 
+1//3 0 0
+-5//48 15//16 0
+]
+b = [1//10, 1//2, 2//5]
+
+
+#up to order 3
+A = [  
+0 0 
+2//3 0 
+]
+b = [1//4, 3//4]
+
+#RK method
+A = [0//1  0//1  0//1  0//1
+ 1//2  0//1  0//1  0//1
+ 0//1  1//2  0//1  0//1
+ 0//1  0//1  1//1  0//1]
+ b =  [1//6, 1//3, 1//3, 1//6]
+
+#this function checks whether a method is energy Preserving for a given order s
+function EnergyPreserving(A,b,s)
+    rka = RungeKuttaMethod(A, b)
+#generate bseries 
+    series_a = modified_equation(bseries(rka, s))
+    coefficients = collect(values(series_a))
+    atrees = collect(keys(series_a))
+# Create an empty vector to store the converted trees
+    trees = Vector{Vector{Int}}(undef, length(series_a))
+# Convert the trees and store them in the trees vector
+    for i in 1:length(series_a)
+        #comment: check Rooted tree object
+        generate_arrays_from_rootedtree = root_converter(atrees[i])
+        if isempty(generate_arrays_from_rootedtree)
+            trees[i] = Int[]
+        else
+            trees[i] = generate_arrays_from_rootedtree
+        end
+    end
+    coefficients = symfact_normalization(coefficients,trees)  
+    signal = IsEnergyPreserving(trees,coefficients)  
+    if signal == false
+        println("Condition Not Satisfied")
+    else
+        println("Condition Satisfied")
+    end
+end
+
+function get_t_arrays(a)
+    t_dict = Dict{Int, Array}()
+    k = a[end]
+    for j in 1:k-1
+        last_j_occurrence = findlast(x -> x == j, a)
+        last_jplus1_occurrence = findlast(x -> x == j+1, a)
+        if isnothing(last_j_occurrence) || isnothing(last_jplus1_occurrence)
+            t_dict[j] = []
+        else
+            t_dict[j] = a[last_j_occurrence+1:last_jplus1_occurrence-1]
+        end
+    end
+    return t_dict
+end
+
+function modify_t_sub(a)
+    t_dict = get_t_arrays(a)
+    m = num_leafs(a)
+    modified_t_dict = Dict{Int, Vector{Int}}()
+    mid_tree = (m+1)/2
+    for j in 1:m
+        if m % 2 == 1 && j == mid_tree
+            modified_t_dict[j] = t_dict[j]
+        elseif m % 2 == 1 && j < mid_tree
+            new_arr = [n+m-2j+1 for n in t_dict[j]]
+            modified_t_dict[j] = new_arr
+        elseif m % 2 == 1 && j > mid_tree
+            new_arr = [m + n - 2*(j) + 1 for n in t_dict[j]]
+            modified_t_dict[j] = new_arr
+        elseif m % 2 == 0
+            if j <= m/2
+                new_arr = [n+m-2j+1 for n in t_dict[j]]
+                modified_t_dict[j] = new_arr
+            elseif j > m/2
+                new_arr = [m + n - 2*(j) + 1  for n in t_dict[j]]
+                modified_t_dict[j] = new_arr
+            end
+        end
+    end
+    return modified_t_dict
+end
+
+function bush_detector(tree)
+    bush = false
+    if tree != [1]
+        l = length(tree)
+        bush = true
+        if l > 1
+            for i in 2:l
+                if tree[i] != 2
+                    bush = false
+                end
+            end
+        end
+    end
+    return bush
+end
+#this function eliminates repeated subarrays
+function eliminate_repeated_subarrays(M)
+    unique_M = []
+    for arr in M
+        if !(arr in unique_M)
+            push!(unique_M, arr)
+        end
+    end
+    return unique_M
+end
+
+
+#this function swaps the trees
+function permuta(a::Vector{Int})
+    diccionario = modify_t_sub(a)
+    m = num_leafs(a)
+    ad_dict = Dict{Int, Vector{Int}}()
+    for j in 1:m
+        ad_dict[j] = diccionario[m-j+1]
+    end
+    return ad_dict
+end
+
+
+
+function adjoint(a::Vector{Int})
+    ad_dict = permuta(a)
+    m = num_leafs(a)
+    adjunto = collect(1:m+1)
+    for j in 1:m
+        last_j_occurrence = findlast(x -> x == j, adjunto)
+        last_jplus1_occurrence = findlast(x -> x == j+1, adjunto)
+        adjunto = vcat(adjunto[1:last_j_occurrence], ad_dict[j], adjunto[last_jplus1_occurrence:end])
+    end
+    adjunto = canonicalarray(adjunto)
+    #println("Adjoint Array: ", adjunto)
+    return adjunto
+end
+
+function root_converter(t::RootedTree{Int64, Vector{Int64}})
+    str = string(t) 
+    arr_str = match(r"\[(.*?)\]", str).captures[1]
+    if isempty(arr_str) || all(iswhitespace, arr_str)
+        return Int[]
+    end
+    arr = parse.(Int, split(arr_str, ","))
+    return arr
+end
+#this function reorders the tree in the same way as the Bseries output does
+function canonicalarray(tree)
+    t = rootedtree(tree)
+    trees = root_converter(t)
+    return trees
+end
+
+#this functions receives
+function indexator(trees,arbol)
+    variable = 0
+    for i in 1:length(trees)
+        if trees[i] == arbol
+            variable = i
+        end
+    end
+    return variable
+end
+
+iswhitespace(c::Char) = isspace(c)
+
+function BMinus(array)
+    m = 1
+    for i in array
+        if i == 2
+            m = m + 1
+        end
+    end
+    contador = 0
+    #aqui almacenamos los 
+    auxiliar = Array{Int64}(undef, m)
+    longitud = length(array)
+    auxiliar[m] = longitud
+    #numero de subtrees
+    n = m-1
+    subtrees = Vector{Any}(undef, n)
+    for i in 2:longitud
+        if array[i] == 2
+            contador += 1
+            auxiliar[contador] = i
+            if contador == n
+                break
+            end
+        end
+    end
+    for i in 1:n
+        subtrees[i] = array[auxiliar[i]:(auxiliar[i+1]-1)]
+    end
+    push!(subtrees[n],array[longitud])
+    return subtrees
+end
+
+#println(BMinus(array))
+
+
+#this function returns a Hamiltonian H matriz of nxn
+
+
+              
+#funcion generadora de arboles simetricos
+function symmetries(tree)
+    subtrees = BMinus(tree)
+    #first order symmetries
+    perms = collect(permutations(subtrees))
+    l = length(perms)
+    for i in 1:l
+        perms[i] = reduce(vcat, perms[i])
+    end
+    result = add_one_to_left(perms)
+    return(result)
+end
+
+#generates all equivalent trees
+function equivalent_trees(array)
+    tree = rootedtree(array)
+    l = length(array)
+    superarray = get_permutations(array)
+    lperm = length(superarray)
+    for i in reverse(1:lperm) 
+        if superarray[i][1] != 1
+            splice!(superarray, i)
+        end
+    end
+    lperm = length(superarray)
+    flag = false
+    for i in reverse(1:lperm)
+        for j in 1:(l-1)
+            if (superarray[i][j+1] - superarray[i][j]) > 1
+                flag = true
+            end
+        end
+        if flag == true
+            splice!(superarray, i)
+        end
+        flag = false
+    end
+    lperm = length(superarray)
+    for i in reverse(1:lperm)
+        arbol = rootedtree(superarray[i])
+        bandera = arbol == tree
+        if bandera == false
+            splice!(superarray, i)
+        end
+    end
+    bandera = false
+    lperm = length(superarray)
+    for i in reverse(1:lperm)
+        for j in 1:(i-1)
+            if superarray[i] == superarray[j]
+                bandera = true
+            end
+        end 
+        if bandera == true
+            splice!(superarray, i)
+        end
+        bandera = false
+    end
+
+    return superarray
+end
+
+function get_permutations(arr)
+    perms = collect(permutations(arr))
+    return perms
+end
+
+function num_leafs(a)
+    k = a[end]
+    return k - 1
+end
+
+function add_one_to_left(arrays)
+    for i in 1:length(arrays)
+        arrays[i] = vcat([1], arrays[i])
+    end
+    return arrays
+end
+
+#this function multplies the coefficient for its symmetry factor
+function symfact_normalization(coef,thetrees)
+    l = length(coef)
+    for i in 1:l
+        factor = 0
+        factor = RootedTrees.symmetry(RootedTree(thetrees[i]))
+        coef[i] = coef[i]*(1//factor)
+    end
+    return coef
+end
+
+#this function tells up to what order a method is Energy Preserving
+function OrderMethod(A,b)
+    s = 1
+    flag = false
+    rka = RungeKuttaMethod(A, b)
+    while flag == false 
+#generate bseries 
+        series_a = modified_equation(bseries(rka, s))
+        coefficients = collect(values(series_a))
+        atrees = collect(keys(series_a))
+# Create an empty vector to store the converted trees
+        trees = Vector{Vector{Int}}(undef, length(series_a))
+# Convert the trees and store them in the trees vector
+        for i in 1:length(series_a)
+            array_from_RTree = root_converter(atrees[i])
+            if isempty(array_from_RTree)
+                trees[i] = Int[]
+            else
+                trees[i] = array_from_RTree
+            end
+        end
+        coefficients = symfact_normalization(coefficients,trees)  
+        #println(trees)
+        #println(coefficients)
+        if IsEnergyPreserving(trees,coefficients) == true
+            flag = true
+        end
+        s = s + 1
+    end
+    println("Energy Preserving for order < ", s-1)
+end
+
+#this function checks if an AVF method is energy Preserving for a given order s
+function EnergyPreservingAVF(s) 
+    series = bseries(s) do t, series
+        if order(t) in (0, 1)
+            return 1 // 1
+        else
+            v = 1 // 1
+            n = 0
+            for subtree in SubtreeIterator(t)
+                v *= series[subtree]
+                n += 1
+            end
+            return v / (n + 1)
+        end
+    end
+#generate bseries 
+    series_a = modified_equation(series)
+    coefficients = collect(values(series_a))
+    atrees = collect(keys(series_a))
+    trees = Vector{Vector{Int}}(undef, length(series_a))
+# Convert the trees and store them in the trees vector
+    for i in 1:length(series_a)
+        array_from_RTree = root_converter(atrees[i])
+       if isempty(array_from_RTree)
+            trees[i] = Int[]
+        else
+            trees[i] = array_from_RTree
+        end
+    end
+    coefficients = symfact_normalization(coefficients,trees)  
+    signal = IsEnergyPreserving(trees,coefficients)  
+    if signal == false
+        println("Condition Not Satisfied")
+    else
+        println("Condition Satisfied")
+    end
+    #series_a
+end
+
+#do not run this function 
+#runs infty
+function OrderAVF() 
+    s = 1
+    flag = false
+    while flag == false 
+        series = bseries(s) do t, series
+            if order(t) in (0, 1)
+                return 1 // 1
+            else
+                v = 1 // 1
+                n = 0
+                for subtree in SubtreeIterator(t)
+                    v *= series[subtree]
+                    n += 1
+                end
+                return v / (n + 1)
+            end
+        end
+#generate bseries 
+        series_a = modified_equation(series)
+        
+        coefficients = collect(values(series_a))
+        atrees = collect(keys(series_a))
+# Create an empty vector to store the converted trees
+        trees = Vector{Vector{Int}}(undef, length(series_a))
+# Convert the trees and store them in the trees vector
+        for i in 1:length(series_a)
+            array_from_RTree = root_converter(atrees[i])
+            if isempty(array_from_RTree)
+                trees[i] = Int[]
+            else
+                trees[i] = array_from_RTree
+            end
+        end
+        coefficients = symfact_normalization(coefficients,trees)  
+        #println(trees)
+        #println(coefficients)
+        if IsEnergyPreserving(trees,coefficients) == false
+            flag = true
+        end
+        s = s + 1
+    end
+    #println("bandera", s)
+    println("Energy Preserving for order < ", s-1)    
+end
+
+
+ 
+function IsEnergyPreserving(trees, coefficients)
+    #obtain the adjoint and check if it exists
+    numero_de_coef = length(trees)
+    MatrizEP = Vector{Int64}[]
+    signal = false
+    #bigflag = false
+    #bandera_array = false
+    for t in trees
+        if !isempty(t) 
+            indice_mayor = indexator(trees,t)
+            #println(indice_mayor)
+            bushflag = bush_detector(t)
+            if bushflag == true
+                if coefficients[indice_mayor] != 0
+                    signal = true
+                end
+                #println("bush")
+            else
+                equiv_set = equivalent_trees(t)
+                bandera_ad = false
+                for arbol in equiv_set
+                #j-th canonical vector
+                    ej = zeros(Int64, numero_de_coef)
+                    ej[indice_mayor] = 1
+                #continue
+                    m = num_leafs(arbol)
+                    t_ad = adjoint(arbol)
+                #revisamos si hay self-adjoint
+                    if (rootedtree(t_ad) == rootedtree(arbol))
+                        bandera_ad = true
+                    end
+                    ek = zeros(Int64, numero_de_coef) 
+                    if t_ad in trees 
+                    #bandera_array = true
+                    #k-th canonical vector
+                        ek[indexator(trees,t_ad)] = 1*(-1)^m
+                    #println(ek)
+                        ej = ej + ek
+                        #println(ej)
+                        push!(MatrizEP, ej)
+                    end
+                end
+            end
+        end
+    end
+    #eliminamos las columnas vacías
+    filter!(x -> any(y -> y != 0, x), MatrizEP)
+    #eliminamos las columnas repetidas
+    MatrizEP = eliminate_repeated_subarrays(MatrizEP)
+    #println(MatrizEP)
+    #println("Longitud de M: ", length(MatrizEP))
+    M = hcat(MatrizEP...)
+    rank_M = rank(M)
+    rank_MV = rank([M coefficients])
+    #println(rank_M == rank_MV)
+    #X = find_vector_X(M,coefficients)
+    #println(X)
+    println("All the trees have been checked:")
+    if signal == true
+        println("Condition Not Satisfied")
+    else
+        println("Condition Satisfied")
+    end
+    return rank_M == rank_MV
+end

From e4d79cca542af8818bbc1c5fe6019f7bdbb0e8f8 Mon Sep 17 00:00:00 2001
From: Sondar74 <gabrielbarriosdeleon@gmail.com>
Date: Mon, 26 Jun 2023 15:25:18 +0300
Subject: [PATCH 02/24] updated 'main' branch and created new branch
 'csrkupdated'. This one should work better than the others.

---
 src/BSeries.jl | 104 ++++++++++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 103 insertions(+), 1 deletion(-)

diff --git a/src/BSeries.jl b/src/BSeries.jl
index fb4f7e2e..134e558c 100644
--- a/src/BSeries.jl
+++ b/src/BSeries.jl
@@ -17,7 +17,7 @@ end
 
 using Latexify: Latexify, LaTeXString
 using Combinatorics: Combinatorics, permutations
-using LinearAlgebra: LinearAlgebra, rank
+using LinearAlgebra: LinearAlgebra, rank, dot
 using SparseArrays: SparseArrays, sparse
 
 @reexport using Polynomials: Polynomials, Polynomial
@@ -40,6 +40,8 @@ export renormalize!
 
 export is_energy_preserving, energy_preserving_order
 
+export elementary_differentials_csrk, CSRK
+
 # Types used for traits
 # These traits may decide between different algorithms based on the
 # corresponding complexity etc.
@@ -74,6 +76,20 @@ end
 
 TruncatedBSeries{T, V}() where {T, V} = TruncatedBSeries{T, V}(OrderedDict{T, V}())
 
+"""
+CSRK struct
+"""
+struct ContinuousStageRungeKuttaMethod{MatT <: AbstractMatrix}
+    matrix::MatT
+end
+
+function CSRK(matrix::AbstractMatrix)
+    T = promote_type(eltype(matrix))
+    _M = T.(matrix)
+    return ContinuousStageRungeKuttaMethod(_M)
+end
+
+
 # general interface methods of `AbstractDict` for `TruncatedBSeries`
 @inline Base.iterate(series::TruncatedBSeries) = iterate(series.coef)
 @inline Base.iterate(series::TruncatedBSeries, state) = iterate(series.coef, state)
@@ -612,6 +628,23 @@ function bseries(ros::RosenbrockMethod, order)
     return series
 end
 
+"""
+bseries CSRK
+"""
+function bseries(csrk::ContinuousStageRungeKuttaMethod, order)
+    csrk = csrk.matrix
+    V = Rational{Int64}
+    series = TruncatedBSeries{RootedTree{Int, Vector{Int}}, V}()
+    series[rootedtree(Int[])] = one(Int64)
+    for o in 1:order
+        for t in RootedTreeIterator(o)
+            series[copy(t)] = elementary_differentials_csrk(csrk, t)
+        end
+    end
+
+    return series
+end
+
 # TODO: bseries(ros::RosenbrockMethod)
 # should create a lazy version, optionally a memoized one
 
@@ -2016,4 +2049,73 @@ function equivalent_trees(tree)
     return equivalent_trees_set
 end
 
+
+
+
+# This function generates a polynomial 
+#       A_{t,z} = [t,t^2/2,..., t^s/s]*M*[1, z, ..., z^(s-1)]^T
+# for a given square matrix M of dimmension s and chars 't' and 'z'.  
+function PolynomialA(M,t,z)
+    s = size(M,1)
+    # we need variables to work with
+    variable1 = Sym(t)
+    variable2 = Sym(z)
+    # conjugate the variable 1, provided that this will be the variable
+    # of the left polynomial
+    variable1 = conjugate(variable1)
+    # generate the components of the polynomial with powers of t
+    poli_z = Array{SymPy.Sym}(undef, s)
+    for i in 1:s
+        poli_z[i] = variable2^(i-1)
+    end
+    # generate the components of the polynomial with powers of z
+    poli_t = Array{SymPy.Sym}(undef, s)
+    for i in 1:s
+        poli_t[i] = (1 // i)*(variable1^i)
+    end
+    # multiply matrix times vector
+    result = M * poli_z
+    # use dot product for the two vectors
+    result = dot(poli_t,result)
+    return result
+end
+
+
+"""
+    elementary_differentials_csrk(M,tree)
+
+This function calculates the CSRK elementary differential for a given 
+square matrix 'M' and a given RootedTree.
+
+"""
+function elementary_differentials_csrk(M,rootedtree)
+    # we'll work with the level sequence
+    tree = rootedtree.level_sequence
+    m = maximum(tree)
+    l = length(tree)
+    # create the variables called 'xi' for 1 <= i <= m
+    variables = []
+    for i in 1:m
+        var_name = "x$i"
+        var = Sym(var_name)
+        push!(variables, var)
+    end
+    inverse_counter = l-1
+    # stablish initial integrand, which is the rightmost leaf (last node of the level sequence)
+    if l > 1
+        integrand = integrate(PolynomialA(M,variables[tree[end]-1],variables[tree[end]]),(variables[tree[end]],0,1))
+    else
+        # if the RootedTree is [1] or [], the elementary differential will be 1.
+        return 1
+    end
+    while inverse_counter > 1
+        pseudo_integrand = PolynomialA(M,variables[tree[inverse_counter]-1],variables[tree[inverse_counter]])*integrand
+        integrand = integrate(pseudo_integrand,(variables[tree[inverse_counter]],0,1))
+        inverse_counter -= 1
+    end
+    # multiply for the Basis_Polynomial, i.e. the Polynomial B
+    # return the integral with respect to x1.
+    return integrate(PolynomialA(M,1,variables[1])*integrand,(variables[1],0,1))
+end
+
 end # module

From 823684d64b722eef619b409a18cb0eaa705d3fca Mon Sep 17 00:00:00 2001
From: Sondar74 <gabrielbarriosdeleon@gmail.com>
Date: Mon, 26 Jun 2023 15:28:02 +0300
Subject: [PATCH 03/24] deleted extra file

---
 Energy_Preserving.jl | 500 -------------------------------------------
 1 file changed, 500 deletions(-)
 delete mode 100644 Energy_Preserving.jl

diff --git a/Energy_Preserving.jl b/Energy_Preserving.jl
deleted file mode 100644
index 4bc481c0..00000000
--- a/Energy_Preserving.jl
+++ /dev/null
@@ -1,500 +0,0 @@
-# Load the packages we will use.
-# These must first be installed using: import Pkg; Pkg.add("package_name")
-#import Pkg; Pkg.add("BSeries")
-using BSeries
-using Latexify  # Only needed for some pretty-printing cells below using `latexify`
-import SymPy; sp=SymPy;
-using LightGraphs
-using Combinatorics
-using RootedTrees
-using LinearAlgebra
-
-#Parametros iniciales
-#stisfies for order <5, and not for 5.
-A = [ 0 0 0 
-1//3 0 0
--5//48 15//16 0
-]
-b = [1//10, 1//2, 2//5]
-
-
-#up to order 3
-A = [  
-0 0 
-2//3 0 
-]
-b = [1//4, 3//4]
-
-#RK method
-A = [0//1  0//1  0//1  0//1
- 1//2  0//1  0//1  0//1
- 0//1  1//2  0//1  0//1
- 0//1  0//1  1//1  0//1]
- b =  [1//6, 1//3, 1//3, 1//6]
-
-#this function checks whether a method is energy Preserving for a given order s
-function EnergyPreserving(A,b,s)
-    rka = RungeKuttaMethod(A, b)
-#generate bseries 
-    series_a = modified_equation(bseries(rka, s))
-    coefficients = collect(values(series_a))
-    atrees = collect(keys(series_a))
-# Create an empty vector to store the converted trees
-    trees = Vector{Vector{Int}}(undef, length(series_a))
-# Convert the trees and store them in the trees vector
-    for i in 1:length(series_a)
-        #comment: check Rooted tree object
-        generate_arrays_from_rootedtree = root_converter(atrees[i])
-        if isempty(generate_arrays_from_rootedtree)
-            trees[i] = Int[]
-        else
-            trees[i] = generate_arrays_from_rootedtree
-        end
-    end
-    coefficients = symfact_normalization(coefficients,trees)  
-    signal = IsEnergyPreserving(trees,coefficients)  
-    if signal == false
-        println("Condition Not Satisfied")
-    else
-        println("Condition Satisfied")
-    end
-end
-
-function get_t_arrays(a)
-    t_dict = Dict{Int, Array}()
-    k = a[end]
-    for j in 1:k-1
-        last_j_occurrence = findlast(x -> x == j, a)
-        last_jplus1_occurrence = findlast(x -> x == j+1, a)
-        if isnothing(last_j_occurrence) || isnothing(last_jplus1_occurrence)
-            t_dict[j] = []
-        else
-            t_dict[j] = a[last_j_occurrence+1:last_jplus1_occurrence-1]
-        end
-    end
-    return t_dict
-end
-
-function modify_t_sub(a)
-    t_dict = get_t_arrays(a)
-    m = num_leafs(a)
-    modified_t_dict = Dict{Int, Vector{Int}}()
-    mid_tree = (m+1)/2
-    for j in 1:m
-        if m % 2 == 1 && j == mid_tree
-            modified_t_dict[j] = t_dict[j]
-        elseif m % 2 == 1 && j < mid_tree
-            new_arr = [n+m-2j+1 for n in t_dict[j]]
-            modified_t_dict[j] = new_arr
-        elseif m % 2 == 1 && j > mid_tree
-            new_arr = [m + n - 2*(j) + 1 for n in t_dict[j]]
-            modified_t_dict[j] = new_arr
-        elseif m % 2 == 0
-            if j <= m/2
-                new_arr = [n+m-2j+1 for n in t_dict[j]]
-                modified_t_dict[j] = new_arr
-            elseif j > m/2
-                new_arr = [m + n - 2*(j) + 1  for n in t_dict[j]]
-                modified_t_dict[j] = new_arr
-            end
-        end
-    end
-    return modified_t_dict
-end
-
-function bush_detector(tree)
-    bush = false
-    if tree != [1]
-        l = length(tree)
-        bush = true
-        if l > 1
-            for i in 2:l
-                if tree[i] != 2
-                    bush = false
-                end
-            end
-        end
-    end
-    return bush
-end
-#this function eliminates repeated subarrays
-function eliminate_repeated_subarrays(M)
-    unique_M = []
-    for arr in M
-        if !(arr in unique_M)
-            push!(unique_M, arr)
-        end
-    end
-    return unique_M
-end
-
-
-#this function swaps the trees
-function permuta(a::Vector{Int})
-    diccionario = modify_t_sub(a)
-    m = num_leafs(a)
-    ad_dict = Dict{Int, Vector{Int}}()
-    for j in 1:m
-        ad_dict[j] = diccionario[m-j+1]
-    end
-    return ad_dict
-end
-
-
-
-function adjoint(a::Vector{Int})
-    ad_dict = permuta(a)
-    m = num_leafs(a)
-    adjunto = collect(1:m+1)
-    for j in 1:m
-        last_j_occurrence = findlast(x -> x == j, adjunto)
-        last_jplus1_occurrence = findlast(x -> x == j+1, adjunto)
-        adjunto = vcat(adjunto[1:last_j_occurrence], ad_dict[j], adjunto[last_jplus1_occurrence:end])
-    end
-    adjunto = canonicalarray(adjunto)
-    #println("Adjoint Array: ", adjunto)
-    return adjunto
-end
-
-function root_converter(t::RootedTree{Int64, Vector{Int64}})
-    str = string(t) 
-    arr_str = match(r"\[(.*?)\]", str).captures[1]
-    if isempty(arr_str) || all(iswhitespace, arr_str)
-        return Int[]
-    end
-    arr = parse.(Int, split(arr_str, ","))
-    return arr
-end
-#this function reorders the tree in the same way as the Bseries output does
-function canonicalarray(tree)
-    t = rootedtree(tree)
-    trees = root_converter(t)
-    return trees
-end
-
-#this functions receives
-function indexator(trees,arbol)
-    variable = 0
-    for i in 1:length(trees)
-        if trees[i] == arbol
-            variable = i
-        end
-    end
-    return variable
-end
-
-iswhitespace(c::Char) = isspace(c)
-
-function BMinus(array)
-    m = 1
-    for i in array
-        if i == 2
-            m = m + 1
-        end
-    end
-    contador = 0
-    #aqui almacenamos los 
-    auxiliar = Array{Int64}(undef, m)
-    longitud = length(array)
-    auxiliar[m] = longitud
-    #numero de subtrees
-    n = m-1
-    subtrees = Vector{Any}(undef, n)
-    for i in 2:longitud
-        if array[i] == 2
-            contador += 1
-            auxiliar[contador] = i
-            if contador == n
-                break
-            end
-        end
-    end
-    for i in 1:n
-        subtrees[i] = array[auxiliar[i]:(auxiliar[i+1]-1)]
-    end
-    push!(subtrees[n],array[longitud])
-    return subtrees
-end
-
-#println(BMinus(array))
-
-
-#this function returns a Hamiltonian H matriz of nxn
-
-
-              
-#funcion generadora de arboles simetricos
-function symmetries(tree)
-    subtrees = BMinus(tree)
-    #first order symmetries
-    perms = collect(permutations(subtrees))
-    l = length(perms)
-    for i in 1:l
-        perms[i] = reduce(vcat, perms[i])
-    end
-    result = add_one_to_left(perms)
-    return(result)
-end
-
-#generates all equivalent trees
-function equivalent_trees(array)
-    tree = rootedtree(array)
-    l = length(array)
-    superarray = get_permutations(array)
-    lperm = length(superarray)
-    for i in reverse(1:lperm) 
-        if superarray[i][1] != 1
-            splice!(superarray, i)
-        end
-    end
-    lperm = length(superarray)
-    flag = false
-    for i in reverse(1:lperm)
-        for j in 1:(l-1)
-            if (superarray[i][j+1] - superarray[i][j]) > 1
-                flag = true
-            end
-        end
-        if flag == true
-            splice!(superarray, i)
-        end
-        flag = false
-    end
-    lperm = length(superarray)
-    for i in reverse(1:lperm)
-        arbol = rootedtree(superarray[i])
-        bandera = arbol == tree
-        if bandera == false
-            splice!(superarray, i)
-        end
-    end
-    bandera = false
-    lperm = length(superarray)
-    for i in reverse(1:lperm)
-        for j in 1:(i-1)
-            if superarray[i] == superarray[j]
-                bandera = true
-            end
-        end 
-        if bandera == true
-            splice!(superarray, i)
-        end
-        bandera = false
-    end
-
-    return superarray
-end
-
-function get_permutations(arr)
-    perms = collect(permutations(arr))
-    return perms
-end
-
-function num_leafs(a)
-    k = a[end]
-    return k - 1
-end
-
-function add_one_to_left(arrays)
-    for i in 1:length(arrays)
-        arrays[i] = vcat([1], arrays[i])
-    end
-    return arrays
-end
-
-#this function multplies the coefficient for its symmetry factor
-function symfact_normalization(coef,thetrees)
-    l = length(coef)
-    for i in 1:l
-        factor = 0
-        factor = RootedTrees.symmetry(RootedTree(thetrees[i]))
-        coef[i] = coef[i]*(1//factor)
-    end
-    return coef
-end
-
-#this function tells up to what order a method is Energy Preserving
-function OrderMethod(A,b)
-    s = 1
-    flag = false
-    rka = RungeKuttaMethod(A, b)
-    while flag == false 
-#generate bseries 
-        series_a = modified_equation(bseries(rka, s))
-        coefficients = collect(values(series_a))
-        atrees = collect(keys(series_a))
-# Create an empty vector to store the converted trees
-        trees = Vector{Vector{Int}}(undef, length(series_a))
-# Convert the trees and store them in the trees vector
-        for i in 1:length(series_a)
-            array_from_RTree = root_converter(atrees[i])
-            if isempty(array_from_RTree)
-                trees[i] = Int[]
-            else
-                trees[i] = array_from_RTree
-            end
-        end
-        coefficients = symfact_normalization(coefficients,trees)  
-        #println(trees)
-        #println(coefficients)
-        if IsEnergyPreserving(trees,coefficients) == true
-            flag = true
-        end
-        s = s + 1
-    end
-    println("Energy Preserving for order < ", s-1)
-end
-
-#this function checks if an AVF method is energy Preserving for a given order s
-function EnergyPreservingAVF(s) 
-    series = bseries(s) do t, series
-        if order(t) in (0, 1)
-            return 1 // 1
-        else
-            v = 1 // 1
-            n = 0
-            for subtree in SubtreeIterator(t)
-                v *= series[subtree]
-                n += 1
-            end
-            return v / (n + 1)
-        end
-    end
-#generate bseries 
-    series_a = modified_equation(series)
-    coefficients = collect(values(series_a))
-    atrees = collect(keys(series_a))
-    trees = Vector{Vector{Int}}(undef, length(series_a))
-# Convert the trees and store them in the trees vector
-    for i in 1:length(series_a)
-        array_from_RTree = root_converter(atrees[i])
-       if isempty(array_from_RTree)
-            trees[i] = Int[]
-        else
-            trees[i] = array_from_RTree
-        end
-    end
-    coefficients = symfact_normalization(coefficients,trees)  
-    signal = IsEnergyPreserving(trees,coefficients)  
-    if signal == false
-        println("Condition Not Satisfied")
-    else
-        println("Condition Satisfied")
-    end
-    #series_a
-end
-
-#do not run this function 
-#runs infty
-function OrderAVF() 
-    s = 1
-    flag = false
-    while flag == false 
-        series = bseries(s) do t, series
-            if order(t) in (0, 1)
-                return 1 // 1
-            else
-                v = 1 // 1
-                n = 0
-                for subtree in SubtreeIterator(t)
-                    v *= series[subtree]
-                    n += 1
-                end
-                return v / (n + 1)
-            end
-        end
-#generate bseries 
-        series_a = modified_equation(series)
-        
-        coefficients = collect(values(series_a))
-        atrees = collect(keys(series_a))
-# Create an empty vector to store the converted trees
-        trees = Vector{Vector{Int}}(undef, length(series_a))
-# Convert the trees and store them in the trees vector
-        for i in 1:length(series_a)
-            array_from_RTree = root_converter(atrees[i])
-            if isempty(array_from_RTree)
-                trees[i] = Int[]
-            else
-                trees[i] = array_from_RTree
-            end
-        end
-        coefficients = symfact_normalization(coefficients,trees)  
-        #println(trees)
-        #println(coefficients)
-        if IsEnergyPreserving(trees,coefficients) == false
-            flag = true
-        end
-        s = s + 1
-    end
-    #println("bandera", s)
-    println("Energy Preserving for order < ", s-1)    
-end
-
-
- 
-function IsEnergyPreserving(trees, coefficients)
-    #obtain the adjoint and check if it exists
-    numero_de_coef = length(trees)
-    MatrizEP = Vector{Int64}[]
-    signal = false
-    #bigflag = false
-    #bandera_array = false
-    for t in trees
-        if !isempty(t) 
-            indice_mayor = indexator(trees,t)
-            #println(indice_mayor)
-            bushflag = bush_detector(t)
-            if bushflag == true
-                if coefficients[indice_mayor] != 0
-                    signal = true
-                end
-                #println("bush")
-            else
-                equiv_set = equivalent_trees(t)
-                bandera_ad = false
-                for arbol in equiv_set
-                #j-th canonical vector
-                    ej = zeros(Int64, numero_de_coef)
-                    ej[indice_mayor] = 1
-                #continue
-                    m = num_leafs(arbol)
-                    t_ad = adjoint(arbol)
-                #revisamos si hay self-adjoint
-                    if (rootedtree(t_ad) == rootedtree(arbol))
-                        bandera_ad = true
-                    end
-                    ek = zeros(Int64, numero_de_coef) 
-                    if t_ad in trees 
-                    #bandera_array = true
-                    #k-th canonical vector
-                        ek[indexator(trees,t_ad)] = 1*(-1)^m
-                    #println(ek)
-                        ej = ej + ek
-                        #println(ej)
-                        push!(MatrizEP, ej)
-                    end
-                end
-            end
-        end
-    end
-    #eliminamos las columnas vacías
-    filter!(x -> any(y -> y != 0, x), MatrizEP)
-    #eliminamos las columnas repetidas
-    MatrizEP = eliminate_repeated_subarrays(MatrizEP)
-    #println(MatrizEP)
-    #println("Longitud de M: ", length(MatrizEP))
-    M = hcat(MatrizEP...)
-    rank_M = rank(M)
-    rank_MV = rank([M coefficients])
-    #println(rank_M == rank_MV)
-    #X = find_vector_X(M,coefficients)
-    #println(X)
-    println("All the trees have been checked:")
-    if signal == true
-        println("Condition Not Satisfied")
-    else
-        println("Condition Satisfied")
-    end
-    return rank_M == rank_MV
-end

From 0df736b3707e55c5eb8373f8585f94f9c2f23119 Mon Sep 17 00:00:00 2001
From: Sondar74 <gabrielbarriosdeleon@gmail.com>
Date: Mon, 26 Jun 2023 15:30:09 +0300
Subject: [PATCH 04/24] fixing 'Spell Check'

---
 src/BSeries.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/BSeries.jl b/src/BSeries.jl
index 134e558c..e750c1a2 100644
--- a/src/BSeries.jl
+++ b/src/BSeries.jl
@@ -2054,7 +2054,7 @@ end
 
 # This function generates a polynomial 
 #       A_{t,z} = [t,t^2/2,..., t^s/s]*M*[1, z, ..., z^(s-1)]^T
-# for a given square matrix M of dimmension s and chars 't' and 'z'.  
+# for a given square matrix M of dimension s and chars 't' and 'z'.  
 function PolynomialA(M,t,z)
     s = size(M,1)
     # we need variables to work with

From 09d4629bbef19aa720fa823698fcc8bd303d34af Mon Sep 17 00:00:00 2001
From: Sondar74 <gabrielbarriosdeleon@gmail.com>
Date: Tue, 27 Jun 2023 11:15:55 +0300
Subject: [PATCH 05/24] added docstrings and removed
 elementary_differentials_csrk from export

---
 src/BSeries.jl | 80 ++++++++++++++++++++++++++++++++++++++------------
 1 file changed, 62 insertions(+), 18 deletions(-)

diff --git a/src/BSeries.jl b/src/BSeries.jl
index e750c1a2..b7d6df7e 100644
--- a/src/BSeries.jl
+++ b/src/BSeries.jl
@@ -40,7 +40,7 @@ export renormalize!
 
 export is_energy_preserving, energy_preserving_order
 
-export elementary_differentials_csrk, CSRK
+export CSRK
 
 # Types used for traits
 # These traits may decide between different algorithms based on the
@@ -76,19 +76,6 @@ end
 
 TruncatedBSeries{T, V}() where {T, V} = TruncatedBSeries{T, V}(OrderedDict{T, V}())
 
-"""
-CSRK struct
-"""
-struct ContinuousStageRungeKuttaMethod{MatT <: AbstractMatrix}
-    matrix::MatT
-end
-
-function CSRK(matrix::AbstractMatrix)
-    T = promote_type(eltype(matrix))
-    _M = T.(matrix)
-    return ContinuousStageRungeKuttaMethod(_M)
-end
-
 
 # general interface methods of `AbstractDict` for `TruncatedBSeries`
 @inline Base.iterate(series::TruncatedBSeries) = iterate(series.coef)
@@ -629,7 +616,61 @@ function bseries(ros::RosenbrockMethod, order)
 end
 
 """
-bseries CSRK
+    ContinuousStageRungeKuttaMethod
+
+A struct that describes a CSRK method. This kind of 'struct' should be constructed 
+via [`CSRK`] 
+        'csrk = CSRK(M)' 
+in order to later call the 'bseries' function. 
+"""
+struct ContinuousStageRungeKuttaMethod{MatT <: AbstractMatrix}
+    matrix::MatT
+end
+
+function CSRK(matrix::AbstractMatrix)
+    T = promote_type(eltype(matrix))
+    _M = T.(matrix)
+    return ContinuousStageRungeKuttaMethod(_M)
+end
+
+"""
+        bseries(csrk::ContinuousStageRungeKuttaMethod, order)
+
+Return a truncated B-series up to the specified `order` with coefficients
+determined by the square matrix located in `csrk`. The coefficients for every 
+RootedTree are obtained via the 'elementary_differentials_csrk' function,
+which is calculated according to Miyatake & Butcher (2015). [See @ref csrk]
+
+# Example:
+The energy-preserving 4x4 matrix given by Miyatake & Butcher (2015) is
+```
+M = [-6//5   72//5  -36//1  24//1;
+        72//5  -144//5 -48//1  72//1;
+        -36//1 -48//1   720//1 -720//1;
+        24//1   72//1  -720//1  720//1]
+```
+
+Then, we calculate the bseries with the following code:
+
+```
+csrk = CSRK(M)
+series = bseries(csrk, 4)
+TruncatedBSeries{RootedTree{Int64, Vector{Int64}}, Rational{Int64}} with 9 entries:
+  RootedTree{Int64}: Int64[]      => 1//1
+  RootedTree{Int64}: [1]          => 1//1
+  RootedTree{Int64}: [1, 2]       => 1//2
+  RootedTree{Int64}: [1, 2, 3]    => 6004799503160661//36028797018963968
+  RootedTree{Int64}: [1, 2, 2]    => 6004799503160661//18014398509481984
+  RootedTree{Int64}: [1, 2, 3, 4] => 6004799503160661//144115188075855872
+  RootedTree{Int64}: [1, 2, 3, 3] => 6004799503160661//72057594037927936
+  RootedTree{Int64}: [1, 2, 3, 2] => 1//8
+  RootedTree{Int64}: [1, 2, 2, 2] => 1//4
+
+```
+# References
+Butcher, John & Miyatake, Yuto. (2015). A Characterization of Energy-Preserving
+Methods and the Construction of Parallel Integrators for Hamiltonian Systems. 
+SIAM Journal on Numerical Analysis. 54. 10.1137/15M1020861. 
 """
 function bseries(csrk::ContinuousStageRungeKuttaMethod, order)
     csrk = csrk.matrix
@@ -2050,8 +2091,6 @@ function equivalent_trees(tree)
 end
 
 
-
-
 # This function generates a polynomial 
 #       A_{t,z} = [t,t^2/2,..., t^s/s]*M*[1, z, ..., z^(s-1)]^T
 # for a given square matrix M of dimension s and chars 't' and 'z'.  
@@ -2085,7 +2124,12 @@ end
     elementary_differentials_csrk(M,tree)
 
 This function calculates the CSRK elementary differential for a given 
-square matrix 'M' and a given RootedTree.
+square matrix 'M' and a given RootedTree according to [@ref].
+
+# References
+Butcher, John & Miyatake, Yuto. (2015). A Characterization of Energy-Preserving
+Methods and the Construction of Parallel Integrators for Hamiltonian Systems. 
+SIAM Journal on Numerical Analysis. 54. 10.1137/15M1020861. 
 
 """
 function elementary_differentials_csrk(M,rootedtree)

From 335265f2674191cb7bd6eeb82d9bb99897b92f04 Mon Sep 17 00:00:00 2001
From: Sondar74 <gabrielbarriosdeleon@gmail.com>
Date: Tue, 27 Jun 2023 12:06:43 +0300
Subject: [PATCH 06/24] modified the comments in the functions 'PolynomialA'
 and 'elementary_differentials_csrk'

---
 src/BSeries.jl | 26 ++++++++++++++++++--------
 1 file changed, 18 insertions(+), 8 deletions(-)

diff --git a/src/BSeries.jl b/src/BSeries.jl
index b7d6df7e..6c55adfe 100644
--- a/src/BSeries.jl
+++ b/src/BSeries.jl
@@ -2095,12 +2095,13 @@ end
 #       A_{t,z} = [t,t^2/2,..., t^s/s]*M*[1, z, ..., z^(s-1)]^T
 # for a given square matrix M of dimension s and chars 't' and 'z'.  
 function PolynomialA(M,t,z)
+    # get the dimension of the matrix
     s = size(M,1)
-    # we need variables to work with
+    # we need symbolic variables to work with
     variable1 = Sym(t)
     variable2 = Sym(z)
-    # conjugate the variable 1, provided that this will be the variable
-    # of the left polynomial
+    # conjugate the variable 1, since  this will be the variable of the left polynomial
+    # and the function 'dot' assumes it to be conjugated
     variable1 = conjugate(variable1)
     # generate the components of the polynomial with powers of t
     poli_z = Array{SymPy.Sym}(undef, s)
@@ -2115,8 +2116,7 @@ function PolynomialA(M,t,z)
     # multiply matrix times vector
     result = M * poli_z
     # use dot product for the two vectors
-    result = dot(poli_t,result)
-    return result
+    return dot(poli_t,result)
 end
 
 
@@ -2130,20 +2130,24 @@ square matrix 'M' and a given RootedTree according to [@ref].
 Butcher, John & Miyatake, Yuto. (2015). A Characterization of Energy-Preserving
 Methods and the Construction of Parallel Integrators for Hamiltonian Systems. 
 SIAM Journal on Numerical Analysis. 54. 10.1137/15M1020861. 
+(https://www.researchgate.net/publication/276211444_A_Characterization_of_Energy-Preserving_Methods_and_the_Construction_of_Parallel_Integrators_for_Hamiltonian_Systems)
 
 """
 function elementary_differentials_csrk(M,rootedtree)
-    # we'll work with the level sequence
+    # we extract the level_sequence of 'rootedtree'
     tree = rootedtree.level_sequence
     m = maximum(tree)
     l = length(tree)
-    # create the variables called 'xi' for 1 <= i <= m
+    # Since we will integrate with symbolic variables for 
+    # every node in the tree, we create the variables called 
+    # 'xi' for 1 <= i <= m, because m is the most distant node from the root. 
     variables = []
     for i in 1:m
         var_name = "x$i"
         var = Sym(var_name)
         push!(variables, var)
     end
+    # we will calculate an integral for every node in the level_sequence from right to left
     inverse_counter = l-1
     # stablish initial integrand, which is the rightmost leaf (last node of the level sequence)
     if l > 1
@@ -2152,12 +2156,18 @@ function elementary_differentials_csrk(M,rootedtree)
         # if the RootedTree is [1] or [], the elementary differential will be 1.
         return 1
     end
+    # Start a cycle for integrating 
     while inverse_counter > 1
+        # we define the pseudo_integrand as the product between the last integral and the new polynomial (since the polynomials are
+        # multiplying each others inside the biggest integral). For a node 'i', this new Polynomial is computed for the variables
+        # 'xi' and 'x(i-1)'. 
         pseudo_integrand = PolynomialA(M,variables[tree[inverse_counter]-1],variables[tree[inverse_counter]])*integrand
+        # integrate this new pseudo_integrand with respect to the variable 'xi'
         integrand = integrate(pseudo_integrand,(variables[tree[inverse_counter]],0,1))
         inverse_counter -= 1
     end
-    # multiply for the Basis_Polynomial, i.e. the Polynomial B
+    # Once we have covered every node except for the base, multiply for the Basis_Polynomial, i.e. the Polynomial B
+    # defined by B_{x1} = A_{1, x1}.
     # return the integral with respect to x1.
     return integrate(PolynomialA(M,1,variables[1])*integrand,(variables[1],0,1))
 end

From 8c4e0b171b4b1e19a6cddab3a8dcaf2907cd0fe9 Mon Sep 17 00:00:00 2001
From: Sondar74 <gabrielbarriosdeleon@gmail.com>
Date: Tue, 27 Jun 2023 13:29:08 +0300
Subject: [PATCH 07/24] Modified the function 'bseries' for CSRK

---
 src/BSeries.jl | 13 +++++--------
 1 file changed, 5 insertions(+), 8 deletions(-)

diff --git a/src/BSeries.jl b/src/BSeries.jl
index 6c55adfe..dc9fe172 100644
--- a/src/BSeries.jl
+++ b/src/BSeries.jl
@@ -674,16 +674,13 @@ SIAM Journal on Numerical Analysis. 54. 10.1137/15M1020861.
 """
 function bseries(csrk::ContinuousStageRungeKuttaMethod, order)
     csrk = csrk.matrix
-    V = Rational{Int64}
-    series = TruncatedBSeries{RootedTree{Int, Vector{Int}}, V}()
-    series[rootedtree(Int[])] = one(Int64)
-    for o in 1:order
-        for t in RootedTreeIterator(o)
-            series[copy(t)] = elementary_differentials_csrk(csrk, t)
+    bseries(o) do t, series
+        if order(t) in (0, 1)
+            return one(csrk)
+        else
+            return elementary_differentials_csrk(csrk, t)
         end
     end
-
-    return series
 end
 
 # TODO: bseries(ros::RosenbrockMethod)

From 7ad53252b6e15e10ab620ea5dc1431e994466a9b Mon Sep 17 00:00:00 2001
From: Sondar74 <gabrielbarriosdeleon@gmail.com>
Date: Tue, 27 Jun 2023 15:47:06 +0300
Subject: [PATCH 08/24] exported 'ContinuousStageRungeKuttaMethod'

---
 src/BSeries.jl | 7 +------
 1 file changed, 1 insertion(+), 6 deletions(-)

diff --git a/src/BSeries.jl b/src/BSeries.jl
index dc9fe172..d5d36835 100644
--- a/src/BSeries.jl
+++ b/src/BSeries.jl
@@ -40,7 +40,7 @@ export renormalize!
 
 export is_energy_preserving, energy_preserving_order
 
-export CSRK
+export ContinuousStageRungeKuttaMethod
 
 # Types used for traits
 # These traits may decide between different algorithms based on the
@@ -627,11 +627,6 @@ struct ContinuousStageRungeKuttaMethod{MatT <: AbstractMatrix}
     matrix::MatT
 end
 
-function CSRK(matrix::AbstractMatrix)
-    T = promote_type(eltype(matrix))
-    _M = T.(matrix)
-    return ContinuousStageRungeKuttaMethod(_M)
-end
 
 """
         bseries(csrk::ContinuousStageRungeKuttaMethod, order)

From 5d82c87316191923d295c87e9cf273f206744380 Mon Sep 17 00:00:00 2001
From: Sondar74 <gabrielbarriosdeleon@gmail.com>
Date: Tue, 27 Jun 2023 18:33:08 +0300
Subject: [PATCH 09/24] added one test, and modified the function 'bseries' for
 CSRK

---
 src/BSeries.jl   | 51 ++++++++++++++++++++++++++++++++++++------------
 test/runtests.jl | 37 +++++++++++++++++++++++++++++++++++
 2 files changed, 75 insertions(+), 13 deletions(-)

diff --git a/src/BSeries.jl b/src/BSeries.jl
index d5d36835..ea7fffde 100644
--- a/src/BSeries.jl
+++ b/src/BSeries.jl
@@ -19,6 +19,7 @@ using Latexify: Latexify, LaTeXString
 using Combinatorics: Combinatorics, permutations
 using LinearAlgebra: LinearAlgebra, rank, dot
 using SparseArrays: SparseArrays, sparse
+using SymPy
 
 @reexport using Polynomials: Polynomials, Polynomial
 
@@ -622,6 +623,13 @@ A struct that describes a CSRK method. This kind of 'struct' should be construct
 via [`CSRK`] 
         'csrk = CSRK(M)' 
 in order to later call the 'bseries' function. 
+
+# References
+- Yuto Miyatake and John C. Butcher.
+"A characterization of energy-preserving methods and the construction of
+parallel integrators for Hamiltonian systems."
+SIAM Journal on Numerical Analysis 54, no. 3 (2016): 
+[DOI: 10.1137/15M1020861](https://doi.org/10.1137/15M1020861) 
 """
 struct ContinuousStageRungeKuttaMethod{MatT <: AbstractMatrix}
     matrix::MatT
@@ -629,12 +637,17 @@ end
 
 
 """
-        bseries(csrk::ContinuousStageRungeKuttaMethod, order)
+    bseries(csrk::ContinuousStageRungeKuttaMethod, order)
 
-Return a truncated B-series up to the specified `order` with coefficients
-determined by the square matrix located in `csrk`. The coefficients for every 
-RootedTree are obtained via the 'elementary_differentials_csrk' function,
-which is calculated according to Miyatake & Butcher (2015). [See @ref csrk]
+Compute the B-series of the [`ContinuousStageRungeKuttaMethod`](@ref) `csrk`
+up to the prescribed integer `order` as described by Miyatake & Butcher (2015).
+    
+!!! note "Normalization by elementary differentials"
+The coefficients of the B-series returned by this method need to be
+multiplied by a power of the time step divided by the `symmetry` of the
+rooted tree and multiplied by the corresponding elementary differential
+of the input vector field ``f``.
+See also [`evaluate`](@ref).
 
 # Example:
 The energy-preserving 4x4 matrix given by Miyatake & Butcher (2015) is
@@ -663,19 +676,31 @@ TruncatedBSeries{RootedTree{Int64, Vector{Int64}}, Rational{Int64}} with 9 entri
 
 ```
 # References
-Butcher, John & Miyatake, Yuto. (2015). A Characterization of Energy-Preserving
-Methods and the Construction of Parallel Integrators for Hamiltonian Systems. 
-SIAM Journal on Numerical Analysis. 54. 10.1137/15M1020861. 
+- Yuto Miyatake and John C. Butcher.
+"A characterization of energy-preserving methods and the construction of
+parallel integrators for Hamiltonian systems."
+SIAM Journal on Numerical Analysis 54, no. 3 (2016): 
+[DOI: 10.1137/15M1020861](https://doi.org/10.1137/15M1020861) 
 """
 function bseries(csrk::ContinuousStageRungeKuttaMethod, order)
     csrk = csrk.matrix
-    bseries(o) do t, series
-        if order(t) in (0, 1)
-            return one(csrk)
-        else
-            return elementary_differentials_csrk(csrk, t)
+    V_tmp = eltype(csrk)
+    if V_tmp <: Integer
+        # If people use integer coefficients, they will likely want to have results
+        # as exact as possible. However, general terms are not integers. Thus, we
+        # use rationals instead.
+        V = Rational{V_tmp}
+    else
+        V = V_tmp
+    end
+    series = TruncatedBSeries{RootedTree{Int, Vector{Int}}, V}()
+    series[rootedtree(Int[])] = one(V)
+    for o in 1:order
+        for t in RootedTreeIterator(o)
+            series[copy(t)] = elementary_differentials_csrk(csrk, t)
         end
     end
+    return series
 end
 
 # TODO: bseries(ros::RosenbrockMethod)
diff --git a/test/runtests.jl b/test/runtests.jl
index 3a6291e2..1290dd78 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2058,4 +2058,41 @@ using Aqua: Aqua
             Aqua.test_project_toml_formatting(BSeries)
         end
     end
+
+    @testset "Continuous Stage Runge-Kutta Method" begin
+        @testset "(Inverse) 4x4-Hilbert Matrix" begin
+            # References
+            # - Yuto Miyatake and John C. Butcher.
+            # "A characterization of energy-preserving methods and the construction of
+            # parallel integrators for Hamiltonian systems."
+            # SIAM Journal on Numerical Analysis 54, no. 3 (2016): 
+            # [DOI: 10.1137/15M1020861](https://doi.org/10.1137/15M1020861) 
+
+            # This is the matrix obtained by inverting the 4x4-Hilbert matrix
+            M = [-6//5   72//5  -36//1  24//1;
+            72//5  -144//5 -48//1  72//1;
+            -36//1 -48//1   720//1 -720//1;
+            24//1   72//1  -720//1  720//1]
+            csrk = ContinuousStageRungeKuttaMethod(M)
+            # Generate the bseries up to order 4
+            order = 4
+            series = bseries(csrk, order)
+            l = length(series)
+            # we save the expected coefficients in the following array
+            expected_coefficients = [1//1 ,1//1, 1//2, 6004799503160661//36028797018963968, 6004799503160661//18014398509481984, 6004799503160661//144115188075855872, 6004799503160661//72057594037927936, 1//8, 1//4]
+            # define a Bool 'coef_match' which checks if all the obtained coefficients are the same as the expected ones.
+            coef_match = true
+            # generate every RootedTree up to order 4, and check if its coefficient in 'series' mathes the expected ones
+            counter = 2
+            for o in 1:order
+                for t in RootedTreeIterator(o)
+                    if series[t] != expected_coefficients[counter]
+                        coef_match = false
+                    end
+                    counter += 1
+                end
+            end
+            @test coef_match == true
+        end
+    end
 end # @testset "BSeries"

From c0a724d1d4db5cf472a2945993b1aee2fd48101d Mon Sep 17 00:00:00 2001
From: Sondar74 <gabrielbarriosdeleon@gmail.com>
Date: Tue, 27 Jun 2023 18:41:08 +0300
Subject: [PATCH 10/24] Spell Check

---
 test/runtests.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index 1290dd78..5748ff37 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2082,7 +2082,7 @@ using Aqua: Aqua
             expected_coefficients = [1//1 ,1//1, 1//2, 6004799503160661//36028797018963968, 6004799503160661//18014398509481984, 6004799503160661//144115188075855872, 6004799503160661//72057594037927936, 1//8, 1//4]
             # define a Bool 'coef_match' which checks if all the obtained coefficients are the same as the expected ones.
             coef_match = true
-            # generate every RootedTree up to order 4, and check if its coefficient in 'series' mathes the expected ones
+            # generate every RootedTree up to order 4, and check if its coefficient in 'series' matches the expected ones
             counter = 2
             for o in 1:order
                 for t in RootedTreeIterator(o)

From c45de1c26422b97ee1cb6c1c8c81a6b58360e348 Mon Sep 17 00:00:00 2001
From: Sondar74 <79122393+Sondar74@users.noreply.github.com>
Date: Wed, 28 Jun 2023 10:33:45 +0300
Subject: [PATCH 11/24] Update src/BSeries.jl

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>
---
 src/BSeries.jl | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/BSeries.jl b/src/BSeries.jl
index ea7fffde..11e91de0 100644
--- a/src/BSeries.jl
+++ b/src/BSeries.jl
@@ -643,11 +643,11 @@ Compute the B-series of the [`ContinuousStageRungeKuttaMethod`](@ref) `csrk`
 up to the prescribed integer `order` as described by Miyatake & Butcher (2015).
     
 !!! note "Normalization by elementary differentials"
-The coefficients of the B-series returned by this method need to be
-multiplied by a power of the time step divided by the `symmetry` of the
-rooted tree and multiplied by the corresponding elementary differential
-of the input vector field ``f``.
-See also [`evaluate`](@ref).
+    The coefficients of the B-series returned by this method need to be
+    multiplied by a power of the time step divided by the `symmetry` of the
+    rooted tree and multiplied by the corresponding elementary differential
+    of the input vector field ``f``.
+    See also [`evaluate`](@ref).
 
 # Example:
 The energy-preserving 4x4 matrix given by Miyatake & Butcher (2015) is

From 6616b9ce04d56e57c234063dbdf783bb25ac2795 Mon Sep 17 00:00:00 2001
From: Sondar74 <79122393+Sondar74@users.noreply.github.com>
Date: Wed, 28 Jun 2023 10:34:44 +0300
Subject: [PATCH 12/24] Update src/BSeries.jl

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>
---
 src/BSeries.jl | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/src/BSeries.jl b/src/BSeries.jl
index 11e91de0..5c33c573 100644
--- a/src/BSeries.jl
+++ b/src/BSeries.jl
@@ -625,11 +625,12 @@ via [`CSRK`]
 in order to later call the 'bseries' function. 
 
 # References
+
 - Yuto Miyatake and John C. Butcher.
-"A characterization of energy-preserving methods and the construction of
-parallel integrators for Hamiltonian systems."
-SIAM Journal on Numerical Analysis 54, no. 3 (2016): 
-[DOI: 10.1137/15M1020861](https://doi.org/10.1137/15M1020861) 
+  "A characterization of energy-preserving methods and the construction of
+  parallel integrators for Hamiltonian systems."
+  SIAM Journal on Numerical Analysis 54, no. 3 (2016): 
+  [DOI: 10.1137/15M1020861](https://doi.org/10.1137/15M1020861) 
 """
 struct ContinuousStageRungeKuttaMethod{MatT <: AbstractMatrix}
     matrix::MatT

From 1227e1ae407a211804ca3966dd380cb256029f5a Mon Sep 17 00:00:00 2001
From: Sondar74 <79122393+Sondar74@users.noreply.github.com>
Date: Wed, 28 Jun 2023 10:35:01 +0300
Subject: [PATCH 13/24] Update src/BSeries.jl

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>
---
 src/BSeries.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/BSeries.jl b/src/BSeries.jl
index 5c33c573..e14ee2ef 100644
--- a/src/BSeries.jl
+++ b/src/BSeries.jl
@@ -659,7 +659,7 @@ M = [-6//5   72//5  -36//1  24//1;
         24//1   72//1  -720//1  720//1]
 ```
 
-Then, we calculate the bseries with the following code:
+Then, we calculate the B-series with the following code:
 
 ```
 csrk = CSRK(M)

From 3691f695096a8b42ea08d3ab26ae88a5cbf4aab1 Mon Sep 17 00:00:00 2001
From: Sondar74 <79122393+Sondar74@users.noreply.github.com>
Date: Wed, 28 Jun 2023 10:40:19 +0300
Subject: [PATCH 14/24] Update src/BSeries.jl

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>
---
 src/BSeries.jl | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/BSeries.jl b/src/BSeries.jl
index e14ee2ef..468c360a 100644
--- a/src/BSeries.jl
+++ b/src/BSeries.jl
@@ -678,10 +678,10 @@ TruncatedBSeries{RootedTree{Int64, Vector{Int64}}, Rational{Int64}} with 9 entri
 ```
 # References
 - Yuto Miyatake and John C. Butcher.
-"A characterization of energy-preserving methods and the construction of
-parallel integrators for Hamiltonian systems."
-SIAM Journal on Numerical Analysis 54, no. 3 (2016): 
-[DOI: 10.1137/15M1020861](https://doi.org/10.1137/15M1020861) 
+  "A characterization of energy-preserving methods and the construction of
+  parallel integrators for Hamiltonian systems."
+  SIAM Journal on Numerical Analysis 54, no. 3 (2016): 
+  [DOI: 10.1137/15M1020861](https://doi.org/10.1137/15M1020861) 
 """
 function bseries(csrk::ContinuousStageRungeKuttaMethod, order)
     csrk = csrk.matrix

From 68cecb02cf6de203c735d979fa06c730f33454ef Mon Sep 17 00:00:00 2001
From: Sondar74 <79122393+Sondar74@users.noreply.github.com>
Date: Wed, 28 Jun 2023 10:43:52 +0300
Subject: [PATCH 15/24] Update test/runtests.jl

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>
---
 test/runtests.jl | 16 +---------------
 1 file changed, 1 insertion(+), 15 deletions(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index 5748ff37..a2c421c8 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2077,22 +2077,8 @@ using Aqua: Aqua
             # Generate the bseries up to order 4
             order = 4
             series = bseries(csrk, order)
-            l = length(series)
-            # we save the expected coefficients in the following array
             expected_coefficients = [1//1 ,1//1, 1//2, 6004799503160661//36028797018963968, 6004799503160661//18014398509481984, 6004799503160661//144115188075855872, 6004799503160661//72057594037927936, 1//8, 1//4]
-            # define a Bool 'coef_match' which checks if all the obtained coefficients are the same as the expected ones.
-            coef_match = true
-            # generate every RootedTree up to order 4, and check if its coefficient in 'series' matches the expected ones
-            counter = 2
-            for o in 1:order
-                for t in RootedTreeIterator(o)
-                    if series[t] != expected_coefficients[counter]
-                        coef_match = false
-                    end
-                    counter += 1
-                end
-            end
-            @test coef_match == true
+            @test collect(values(series)) == expected_coefficients
         end
     end
 end # @testset "BSeries"

From 5d3086595d60420e2eee259e7d3964c701742154 Mon Sep 17 00:00:00 2001
From: Sondar74 <gabrielbarriosdeleon@gmail.com>
Date: Wed, 28 Jun 2023 10:45:27 +0300
Subject: [PATCH 16/24] updating docstring

---
 src/BSeries.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/BSeries.jl b/src/BSeries.jl
index ea7fffde..8b579cd5 100644
--- a/src/BSeries.jl
+++ b/src/BSeries.jl
@@ -619,9 +619,9 @@ end
 """
     ContinuousStageRungeKuttaMethod
 
-A struct that describes a CSRK method. This kind of 'struct' should be constructed 
-via [`CSRK`] 
-        'csrk = CSRK(M)' 
+A struct that describes a CSRK method. This kind of `struct` should be constructed 
+via [ContinuousStageRungeKuttaMethod] 
+    `csrk = ContinuousStageRungeKuttaMethod(M)`
 in order to later call the 'bseries' function. 
 
 # References

From 53fd493de4ac441c594f8277bbe1838ec2b957a3 Mon Sep 17 00:00:00 2001
From: Sondar74 <gabrielbarriosdeleon@gmail.com>
Date: Wed, 28 Jun 2023 11:47:51 +0300
Subject: [PATCH 17/24] fixed the values of bseries()

---
 src/BSeries.jl   | 10 +++++-----
 test/runtests.jl |  2 +-
 2 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/src/BSeries.jl b/src/BSeries.jl
index 036db735..6bf1bfb3 100644
--- a/src/BSeries.jl
+++ b/src/BSeries.jl
@@ -668,10 +668,10 @@ TruncatedBSeries{RootedTree{Int64, Vector{Int64}}, Rational{Int64}} with 9 entri
   RootedTree{Int64}: Int64[]      => 1//1
   RootedTree{Int64}: [1]          => 1//1
   RootedTree{Int64}: [1, 2]       => 1//2
-  RootedTree{Int64}: [1, 2, 3]    => 6004799503160661//36028797018963968
-  RootedTree{Int64}: [1, 2, 2]    => 6004799503160661//18014398509481984
-  RootedTree{Int64}: [1, 2, 3, 4] => 6004799503160661//144115188075855872
-  RootedTree{Int64}: [1, 2, 3, 3] => 6004799503160661//72057594037927936
+  RootedTree{Int64}: [1, 2, 3]    => 1//6
+  RootedTree{Int64}: [1, 2, 2]    => 1//3
+  RootedTree{Int64}: [1, 2, 3, 4] => 1//24
+  RootedTree{Int64}: [1, 2, 3, 3] => 1//12
   RootedTree{Int64}: [1, 2, 3, 2] => 1//8
   RootedTree{Int64}: [1, 2, 2, 2] => 1//4
 
@@ -698,7 +698,7 @@ function bseries(csrk::ContinuousStageRungeKuttaMethod, order)
     series[rootedtree(Int[])] = one(V)
     for o in 1:order
         for t in RootedTreeIterator(o)
-            series[copy(t)] = elementary_differentials_csrk(csrk, t)
+            series[copy(t)] = N(elementary_differentials_csrk(csrk, t))
         end
     end
     return series
diff --git a/test/runtests.jl b/test/runtests.jl
index a2c421c8..8351040e 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2077,7 +2077,7 @@ using Aqua: Aqua
             # Generate the bseries up to order 4
             order = 4
             series = bseries(csrk, order)
-            expected_coefficients = [1//1 ,1//1, 1//2, 6004799503160661//36028797018963968, 6004799503160661//18014398509481984, 6004799503160661//144115188075855872, 6004799503160661//72057594037927936, 1//8, 1//4]
+            expected_coefficients = [1//1 ,1//1, 1//2, 1//6, 1//3, 1//24, 1//12, 1//8, 1//4]
             @test collect(values(series)) == expected_coefficients
         end
     end

From d6a23d6513d5b9e66b35393f2449f151a33b5503 Mon Sep 17 00:00:00 2001
From: Sondar74 <gabrielbarriosdeleon@gmail.com>
Date: Wed, 28 Jun 2023 13:34:06 +0300
Subject: [PATCH 18/24] added test for floating-point and symbolic
 coefficientes using SymPy and SymEngine

---
 test/runtests.jl | 37 +++++++++++++++++++++++++++++++++++++
 1 file changed, 37 insertions(+)

diff --git a/test/runtests.jl b/test/runtests.jl
index 8351040e..6972cde7 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2080,5 +2080,42 @@ using Aqua: Aqua
             expected_coefficients = [1//1 ,1//1, 1//2, 1//6, 1//3, 1//24, 1//12, 1//8, 1//4]
             @test collect(values(series)) == expected_coefficients
         end
+
+        @testset "Floating-points coefficients" begin
+            # Define variables
+            a = -0.37069987
+            b = 0.74139974
+            c = 2.51720052
+
+            # Create symbolic matrix
+            Minv = [a 1//2 b 1//6;
+            1//2 1//4 1//6 1//8;
+            b 1//6 c 1//10;
+            1//6 1//8 1//10 1//12]
+
+            # Calculate inverse of the  matrix
+            M = inv(Minv)
+            csrk = ContinuousStageRungeKuttaMethod(M)
+            # Generate the bseries up to order 4
+            order = 4
+            series = bseries(csrk, order)
+            expected_coefficients = [1.0 ,1.0, 0.5000000000000002, 0.16666666666666666, 0.3333333333333336, 0.04166666666666661, 0.0833333333333332, 0.12500000000000003, 0.2500000000000003]
+            @test collect(values(series)) == expected_coefficients
+        end
+        @testset "Symbolic coefficients using SymPy.jl" begin
+            # Define variables
+            import SymPy: symbols
+            x = symbols("x")
+            y = symbols("y")
+            # Create symbolic matrix
+            M = [x y;
+                y x]
+            csrk = ContinuousStageRungeKuttaMethod(M)
+            # Generate the bseries up to order 4
+            order = 3
+            series = bseries(csrk, order)
+            expected_coefficients = [1,1, 25*x^2/32 + 5*x*y/4 + y^2/2,  1105*x^3/2304 + 671*x^2*y/576 + 545*x*y^2/576 + 37*y^3/144, 125*x^3/192 + 25*x^2*y/16 + 5*x*y^2/4 + y^3/3]
+            @test collect(values(series)) == expected_coefficients
+        end
     end
 end # @testset "BSeries"

From c0ba86fe242058e62e1f63eadb7dcf3fe3b54539 Mon Sep 17 00:00:00 2001
From: Sondar74 <gabrielbarriosdeleon@gmail.com>
Date: Wed, 28 Jun 2023 14:16:11 +0300
Subject: [PATCH 19/24] added test for Symbolics.jl: same as SymEngine

---
 test/runtests.jl | 28 ++++++++++++++++++++++++++++
 1 file changed, 28 insertions(+)

diff --git a/test/runtests.jl b/test/runtests.jl
index 6972cde7..f1dbb529 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2117,5 +2117,33 @@ using Aqua: Aqua
             expected_coefficients = [1,1, 25*x^2/32 + 5*x*y/4 + y^2/2,  1105*x^3/2304 + 671*x^2*y/576 + 545*x*y^2/576 + 37*y^3/144, 125*x^3/192 + 25*x^2*y/16 + 5*x*y^2/4 + y^3/3]
             @test collect(values(series)) == expected_coefficients
         end
+        @testset "Symbolic coefficients using SymEngine.jl" begin
+            # Define variables
+            import SymEngine: symbols
+            x,y = SymEngine.symbols("x y")
+            # Create symbolic matrix
+            M = [x y;
+                y x]
+            csrk = ContinuousStageRungeKuttaMethod(M)
+            # Generate the bseries up to order 4
+            order = 3
+            series = bseries(csrk, order)
+            expected_coefficients = [1,1, 25*x^2/32 + 5*x*y/4 + y^2/2,  1105*x^3/2304 + 671*x^2*y/576 + 545*x*y^2/576 + 37*y^3/144, 125*x^3/192 + 25*x^2*y/16 + 5*x*y^2/4 + y^3/3]
+            @test collect(values(series)) == expected_coefficients
+        end
+        @testset "Symbolic coefficients using Symbolics.jl" begin
+            # Define variables
+            import Symbolics: @variables
+            @variables x y
+            # Create symbolic matrix
+            M = [x y;
+                y x]
+            csrk = ContinuousStageRungeKuttaMethod(M)
+            # Generate the bseries up to order 4
+            order = 3
+            series = bseries(csrk, order)
+            expected_coefficients = [1,1, 25*x^2/32 + 5*x*y/4 + y^2/2,  1105*x^3/2304 + 671*x^2*y/576 + 545*x*y^2/576 + 37*y^3/144, 125*x^3/192 + 25*x^2*y/16 + 5*x*y^2/4 + y^3/3]
+            @test collect(values(series)) == expected_coefficients
+        end
     end
 end # @testset "BSeries"

From b5813d263e8d167771d2739c059049446d699ec7 Mon Sep 17 00:00:00 2001
From: Sondar74 <79122393+Sondar74@users.noreply.github.com>
Date: Wed, 28 Jun 2023 17:06:19 +0300
Subject: [PATCH 20/24] Update test/runtests.jl

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>
---
 test/runtests.jl | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index f1dbb529..ef81ae16 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2133,8 +2133,7 @@ using Aqua: Aqua
         end
         @testset "Symbolic coefficients using Symbolics.jl" begin
             # Define variables
-            import Symbolics: @variables
-            @variables x y
+            Symbolics.@variables x y
             # Create symbolic matrix
             M = [x y;
                 y x]

From 691507f5dee87f57cb1dec9a8d739342b6b9d76f Mon Sep 17 00:00:00 2001
From: Sondar74 <79122393+Sondar74@users.noreply.github.com>
Date: Wed, 28 Jun 2023 17:06:30 +0300
Subject: [PATCH 21/24] Update test/runtests.jl

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>
---
 test/runtests.jl | 1 -
 1 file changed, 1 deletion(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index ef81ae16..46baa694 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2119,7 +2119,6 @@ using Aqua: Aqua
         end
         @testset "Symbolic coefficients using SymEngine.jl" begin
             # Define variables
-            import SymEngine: symbols
             x,y = SymEngine.symbols("x y")
             # Create symbolic matrix
             M = [x y;

From 02acff4577dca84ac6cda06fb55ab2d86061599f Mon Sep 17 00:00:00 2001
From: Sondar74 <79122393+Sondar74@users.noreply.github.com>
Date: Wed, 28 Jun 2023 17:12:10 +0300
Subject: [PATCH 22/24] Update test/runtests.jl

Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com>
---
 test/runtests.jl | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index 46baa694..60791a04 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2104,9 +2104,8 @@ using Aqua: Aqua
         end
         @testset "Symbolic coefficients using SymPy.jl" begin
             # Define variables
-            import SymPy: symbols
-            x = symbols("x")
-            y = symbols("y")
+            x = SymPy.symbols("x")
+            y = SymPy.symbols("y")
             # Create symbolic matrix
             M = [x y;
                 y x]

From fd1a23e07b4726c037f9ab2c5072bd502a101093 Mon Sep 17 00:00:00 2001
From: Sondar74 <gabrielbarriosdeleon@gmail.com>
Date: Wed, 28 Jun 2023 17:46:33 +0300
Subject: [PATCH 23/24] added Energy_Preserving test

---
 src/BSeries.jl   |  6 ++++--
 test/runtests.jl | 38 +++++++++++++++++++++++++++++++++-----
 2 files changed, 37 insertions(+), 7 deletions(-)

diff --git a/src/BSeries.jl b/src/BSeries.jl
index 6bf1bfb3..d68e9955 100644
--- a/src/BSeries.jl
+++ b/src/BSeries.jl
@@ -632,10 +632,12 @@ in order to later call the 'bseries' function.
   SIAM Journal on Numerical Analysis 54, no. 3 (2016): 
   [DOI: 10.1137/15M1020861](https://doi.org/10.1137/15M1020861) 
 """
-struct ContinuousStageRungeKuttaMethod{MatT <: AbstractMatrix}
-    matrix::MatT
+struct ContinuousStageRungeKuttaMethod{T, MatT <: AbstractMatrix{T}} <: RootedTrees.AbstractTimeIntegrationMethod
+    A::MatT
 end
 
+Base.eltype(::ContinuousStageRungeKuttaMethod{T}) where {T} = T
+
 
 """
     bseries(csrk::ContinuousStageRungeKuttaMethod, order)
diff --git a/test/runtests.jl b/test/runtests.jl
index f1dbb529..6b3def16 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2083,11 +2083,19 @@ using Aqua: Aqua
 
         @testset "Floating-points coefficients" begin
             # Define variables
+            # Values from current research for CSRK Methods.
+            
+            # References
+            # - Yuto Miyatake and John C. Butcher.
+            # "A characterization of energy-preserving methods and the construction of
+            # parallel integrators for Hamiltonian systems."
+            # SIAM Journal on Numerical Analysis 54, no. 3 (2016): 
+            # [DOI: 10.1137/15M1020861](https://doi.org/10.1137/15M1020861)
             a = -0.37069987
             b = 0.74139974
             c = 2.51720052
 
-            # Create symbolic matrix
+            # Create the 4x4 matrix
             Minv = [a 1//2 b 1//6;
             1//2 1//4 1//6 1//8;
             b 1//6 c 1//10;
@@ -2111,7 +2119,7 @@ using Aqua: Aqua
             M = [x y;
                 y x]
             csrk = ContinuousStageRungeKuttaMethod(M)
-            # Generate the bseries up to order 4
+            # Generate the bseries up to order 3
             order = 3
             series = bseries(csrk, order)
             expected_coefficients = [1,1, 25*x^2/32 + 5*x*y/4 + y^2/2,  1105*x^3/2304 + 671*x^2*y/576 + 545*x*y^2/576 + 37*y^3/144, 125*x^3/192 + 25*x^2*y/16 + 5*x*y^2/4 + y^3/3]
@@ -2125,7 +2133,7 @@ using Aqua: Aqua
             M = [x y;
                 y x]
             csrk = ContinuousStageRungeKuttaMethod(M)
-            # Generate the bseries up to order 4
+            # Generate the bseries up to order 3
             order = 3
             series = bseries(csrk, order)
             expected_coefficients = [1,1, 25*x^2/32 + 5*x*y/4 + y^2/2,  1105*x^3/2304 + 671*x^2*y/576 + 545*x*y^2/576 + 37*y^3/144, 125*x^3/192 + 25*x^2*y/16 + 5*x*y^2/4 + y^3/3]
@@ -2134,16 +2142,36 @@ using Aqua: Aqua
         @testset "Symbolic coefficients using Symbolics.jl" begin
             # Define variables
             import Symbolics: @variables
-            @variables x y
+            (Symbolics.@variables x y)
             # Create symbolic matrix
             M = [x y;
                 y x]
             csrk = ContinuousStageRungeKuttaMethod(M)
-            # Generate the bseries up to order 4
+            # Generate the bseries up to order 3
             order = 3
             series = bseries(csrk, order)
             expected_coefficients = [1,1, 25*x^2/32 + 5*x*y/4 + y^2/2,  1105*x^3/2304 + 671*x^2*y/576 + 545*x*y^2/576 + 37*y^3/144, 125*x^3/192 + 25*x^2*y/16 + 5*x*y^2/4 + y^3/3]
             @test collect(values(series)) == expected_coefficients
         end
+        @testset "Energy-Preserving CSRK" begin
+            # References
+            # - Yuto Miyatake and John C. Butcher.
+            # "A characterization of energy-preserving methods and the construction of
+            # parallel integrators for Hamiltonian systems."
+            # SIAM Journal on Numerical Analysis 54, no. 3 (2016): 
+            # [DOI: 10.1137/15M1020861](https://doi.org/10.1137/15M1020861) 
+
+            # This is the matrix obtained by inverting the 4x4-Hilbert matrix
+            M = [-6//5   72//5  -36//1  24//1;
+            72//5  -144//5 -48//1  72//1;
+            -36//1 -48//1   720//1 -720//1;
+            24//1   72//1  -720//1  720//1]
+            csrk = ContinuousStageRungeKuttaMethod(M)
+            # Generate the bseries up to order 5
+            order = 5
+            series = bseries(csrk, order)
+            #check if it is energy-preserving
+            @test is_energy_preserving(series) == true
+        end
     end
 end # @testset "BSeries"

From c8113e6f1b0b4fcdbb1f064981df7db42beaf533 Mon Sep 17 00:00:00 2001
From: Sondar74 <gabrielbarriosdeleon@gmail.com>
Date: Mon, 10 Jul 2023 16:03:37 +0300
Subject: [PATCH 24/24] fixed the funtion bseries(csrk, order)

---
 src/BSeries.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/BSeries.jl b/src/BSeries.jl
index d68e9955..0872ad67 100644
--- a/src/BSeries.jl
+++ b/src/BSeries.jl
@@ -632,7 +632,7 @@ in order to later call the 'bseries' function.
   SIAM Journal on Numerical Analysis 54, no. 3 (2016): 
   [DOI: 10.1137/15M1020861](https://doi.org/10.1137/15M1020861) 
 """
-struct ContinuousStageRungeKuttaMethod{T, MatT <: AbstractMatrix{T}} <: RootedTrees.AbstractTimeIntegrationMethod
+struct ContinuousStageRungeKuttaMethod{MatT <: AbstractMatrix} <: RootedTrees.AbstractTimeIntegrationMethod
     A::MatT
 end
 
@@ -664,7 +664,7 @@ M = [-6//5   72//5  -36//1  24//1;
 Then, we calculate the B-series with the following code:
 
 ```
-csrk = CSRK(M)
+csrk = ContinuousStageRungeKuttaMethod(M)
 series = bseries(csrk, 4)
 TruncatedBSeries{RootedTree{Int64, Vector{Int64}}, Rational{Int64}} with 9 entries:
   RootedTree{Int64}: Int64[]      => 1//1
@@ -686,7 +686,7 @@ TruncatedBSeries{RootedTree{Int64, Vector{Int64}}, Rational{Int64}} with 9 entri
   [DOI: 10.1137/15M1020861](https://doi.org/10.1137/15M1020861) 
 """
 function bseries(csrk::ContinuousStageRungeKuttaMethod, order)
-    csrk = csrk.matrix
+    csrk = csrk.A
     V_tmp = eltype(csrk)
     if V_tmp <: Integer
         # If people use integer coefficients, they will likely want to have results