From bf98fefc2efca2935aed14ac0e794dac873be53e Mon Sep 17 00:00:00 2001
From: Mikel <mikel.antonana@gmail.com>
Date: Thu, 28 Mar 2024 17:43:07 +0100
Subject: [PATCH] issue 61

---
 ...ti Differencial Equations-checkpoint.ipynb | 236 ++++++++++++++++++
 Examples/Riccati Differencial Equations.ipynb | 236 ++++++++++++++++++
 ODEProblems/InitialNBody15.jl                 |  26 +-
 ODEProblems/InitialNBody16.jl                 |  30 +--
 ODEProblems/InitialNBody9.jl                  |   6 +-
 ODEProblems/InitialNLS.jl                     |   4 +-
 ODEProblems/InitialPleiades.jl                |   4 +-
 src/IRKCoefficients.jl                        |  12 +-
 src/IRKGL16AuxFunctions.jl                    |   4 +
 src/IRKGL16Solver.jl                          |  19 +-
 test/runtests.jl                              |   1 -
 11 files changed, 527 insertions(+), 51 deletions(-)
 create mode 100644 Examples/.ipynb_checkpoints/Riccati Differencial Equations-checkpoint.ipynb
 create mode 100644 Examples/Riccati Differencial Equations.ipynb

diff --git a/Examples/.ipynb_checkpoints/Riccati Differencial Equations-checkpoint.ipynb b/Examples/.ipynb_checkpoints/Riccati Differencial Equations-checkpoint.ipynb
new file mode 100644
index 0000000..5dd6f71
--- /dev/null
+++ b/Examples/.ipynb_checkpoints/Riccati Differencial Equations-checkpoint.ipynb	
@@ -0,0 +1,236 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2fb54ab2",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "using IRKGaussLegendre\n",
+    "using OrdinaryDiffEq"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ce474231",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "function x(t)\n",
+    "    u0 = [1.]\n",
+    "    tspan = (1.,t)\n",
+    "    prob = ODEProblem(Ric1ODE!, u0, tspan)\n",
+    "    sol = solve(prob, Tsit5(); reltol = 1.e-10, abstol = 1.e-10, save_everystep = false)\n",
+    "    return sol(t)\n",
+    "end\n",
+    "\n",
+    "function Ric1ODE!(du, u, pars, t) \n",
+    "    du[1] = -2u[1]-1+u[1]^2\n",
+    "end"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4c16810d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "f(u, p, t) = x(t)[1]^2\n",
+    "\n",
+    "function f1!(du, u, pars, t)\n",
+    "    du[1] = pars(t)\n",
+    "end\n",
+    "\n",
+    "u0 = [0.]; u00 = 0.;\n",
+    "tspan = (0.0, 1.0)\n",
+    "pars = t->x(t)[1]^2\n",
+    "prob = ODEProblem(f, u00, tspan)\n",
+    "prob1 = ODEProblem(f1!, u0, tspan, pars)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "edd522be",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sol1 = solve(prob1, Tsit5(), reltol = 1e-8, abstol = 1e-8); sol1(1)[1]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "899dbbb7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "eltype(prob1.tspan)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "dc97f9b9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sol2 = solve(prob1, IRKGaussLegendre.IRKGL16(maxtrials=4), dt=1., reltol = 1e-8, abstol = 1e-8); sol2(1)[1]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "5ce7b8ef",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# v 1.10.0 :  3.3006156872882713"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "3c71fc4e",
+   "metadata": {},
+   "source": [
+    "## Andreas examples"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b8ccc3cc",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sol2 = solve(prob1, IRKGaussLegendre.IRKGL16(), reltol = 1e-8, abstol = 1e-8); sol2(1)[1]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ab14de53",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sol2 = solve(prob1, IRKGaussLegendre.IRKGL16(), dt = 1e-7,reltol = 1e-8, abstol = 1e-8); sol2(1)[1]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "50dce94e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sol2 = solve(prob1, IRKGaussLegendre.IRKGL16(), dt =1e-6,reltol = 1e-8, abstol = 1e-8); sol2(1)[1]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "fef268bc",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sol2 = solve(prob1, IRKGaussLegendre.IRKGL16(), dt = 1e-5,reltol = 1e-8, abstol = 1e-8); sol2(1)[1]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "86e473b0",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sol2 = solve(prob1, IRKGaussLegendre.IRKGL16(), dt = 1e-4,reltol = 1e-8, abstol = 1e-8); sol2(1)[1]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e15976ba",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sol2 = solve(prob1, IRKGaussLegendre.IRKGL16(), dt = 1e-3,reltol = 1e-8, abstol = 1e-8); sol2(1)[1]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "fe98afa7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sol2 = solve(prob1, IRKGaussLegendre.IRKGL16(), dt = 1e-2,reltol = 1e-8, abstol = 1e-8); sol2(1)[1]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3f241056",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sol2 = solve(prob1, IRKGaussLegendre.IRKGL16(), dt = 1e-1,reltol = 1e-8, abstol = 1e-8); sol2(1)[1]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0d05cf95",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sol2 = solve(prob1, IRKGaussLegendre.IRKGL16(), dt = 1,reltol = 1e-8, abstol = 1e-8); sol2(1)[1]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b4a7cf4a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sol2 = solve(prob1, IRKGaussLegendre.IRKGL16(), dt = 10,reltol = 1e-8, abstol = 1e-8); sol2(1)[1]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "cd3ee271",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sol2 = solve(prob1, IRKGaussLegendre.IRKGL16(), dt = 100,reltol = 1e-8, abstol = 1e-8); sol2(1)[1]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d18c366f",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Julia 1.9.1",
+   "language": "julia",
+   "name": "julia-1.9"
+  },
+  "language_info": {
+   "file_extension": ".jl",
+   "mimetype": "application/julia",
+   "name": "julia",
+   "version": "1.9.1"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/Examples/Riccati Differencial Equations.ipynb b/Examples/Riccati Differencial Equations.ipynb
new file mode 100644
index 0000000..5dd6f71
--- /dev/null
+++ b/Examples/Riccati Differencial Equations.ipynb	
@@ -0,0 +1,236 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2fb54ab2",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "using IRKGaussLegendre\n",
+    "using OrdinaryDiffEq"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ce474231",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "function x(t)\n",
+    "    u0 = [1.]\n",
+    "    tspan = (1.,t)\n",
+    "    prob = ODEProblem(Ric1ODE!, u0, tspan)\n",
+    "    sol = solve(prob, Tsit5(); reltol = 1.e-10, abstol = 1.e-10, save_everystep = false)\n",
+    "    return sol(t)\n",
+    "end\n",
+    "\n",
+    "function Ric1ODE!(du, u, pars, t) \n",
+    "    du[1] = -2u[1]-1+u[1]^2\n",
+    "end"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4c16810d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "f(u, p, t) = x(t)[1]^2\n",
+    "\n",
+    "function f1!(du, u, pars, t)\n",
+    "    du[1] = pars(t)\n",
+    "end\n",
+    "\n",
+    "u0 = [0.]; u00 = 0.;\n",
+    "tspan = (0.0, 1.0)\n",
+    "pars = t->x(t)[1]^2\n",
+    "prob = ODEProblem(f, u00, tspan)\n",
+    "prob1 = ODEProblem(f1!, u0, tspan, pars)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "edd522be",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sol1 = solve(prob1, Tsit5(), reltol = 1e-8, abstol = 1e-8); sol1(1)[1]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "899dbbb7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "eltype(prob1.tspan)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "dc97f9b9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sol2 = solve(prob1, IRKGaussLegendre.IRKGL16(maxtrials=4), dt=1., reltol = 1e-8, abstol = 1e-8); sol2(1)[1]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "5ce7b8ef",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# v 1.10.0 :  3.3006156872882713"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "3c71fc4e",
+   "metadata": {},
+   "source": [
+    "## Andreas examples"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b8ccc3cc",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sol2 = solve(prob1, IRKGaussLegendre.IRKGL16(), reltol = 1e-8, abstol = 1e-8); sol2(1)[1]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ab14de53",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sol2 = solve(prob1, IRKGaussLegendre.IRKGL16(), dt = 1e-7,reltol = 1e-8, abstol = 1e-8); sol2(1)[1]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "50dce94e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sol2 = solve(prob1, IRKGaussLegendre.IRKGL16(), dt =1e-6,reltol = 1e-8, abstol = 1e-8); sol2(1)[1]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "fef268bc",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sol2 = solve(prob1, IRKGaussLegendre.IRKGL16(), dt = 1e-5,reltol = 1e-8, abstol = 1e-8); sol2(1)[1]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "86e473b0",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sol2 = solve(prob1, IRKGaussLegendre.IRKGL16(), dt = 1e-4,reltol = 1e-8, abstol = 1e-8); sol2(1)[1]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e15976ba",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sol2 = solve(prob1, IRKGaussLegendre.IRKGL16(), dt = 1e-3,reltol = 1e-8, abstol = 1e-8); sol2(1)[1]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "fe98afa7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sol2 = solve(prob1, IRKGaussLegendre.IRKGL16(), dt = 1e-2,reltol = 1e-8, abstol = 1e-8); sol2(1)[1]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3f241056",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sol2 = solve(prob1, IRKGaussLegendre.IRKGL16(), dt = 1e-1,reltol = 1e-8, abstol = 1e-8); sol2(1)[1]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0d05cf95",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sol2 = solve(prob1, IRKGaussLegendre.IRKGL16(), dt = 1,reltol = 1e-8, abstol = 1e-8); sol2(1)[1]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b4a7cf4a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sol2 = solve(prob1, IRKGaussLegendre.IRKGL16(), dt = 10,reltol = 1e-8, abstol = 1e-8); sol2(1)[1]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "cd3ee271",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sol2 = solve(prob1, IRKGaussLegendre.IRKGL16(), dt = 100,reltol = 1e-8, abstol = 1e-8); sol2(1)[1]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d18c366f",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Julia 1.9.1",
+   "language": "julia",
+   "name": "julia-1.9"
+  },
+  "language_info": {
+   "file_extension": ".jl",
+   "mimetype": "application/julia",
+   "name": "julia",
+   "version": "1.9.1"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/ODEProblems/InitialNBody15.jl b/ODEProblems/InitialNBody15.jl
index 10163ac..507f033 100644
--- a/ODEProblems/InitialNBody15.jl
+++ b/ODEProblems/InitialNBody15.jl
@@ -22,7 +22,7 @@ function InitialNBody15(T)
         parse(BigFloat, "0.310444819893871300e-13"),
         parse(BigFloat, "0.385475018780881000e-13"),
         parse(BigFloat, "0.213643444257140700e-14"),
-        parse(BigFloat, "0.138862658985619900e-14"),
+        parse(BigFloat, "0.138862658985619900e-14")
     ]
 
     N = length(Gm)
@@ -57,7 +57,7 @@ function InitialNBody15(T)
         parse(BigFloat, "-9.40015796880239402236"),
         parse(BigFloat, "-30.48331376718383944523"),
         parse(BigFloat, "-0.87240555684104999295"),
-        parse(BigFloat, "8.91157617249954997764"),
+        parse(BigFloat, "8.91157617249954997764")
     ]
 
     v = [
@@ -90,7 +90,7 @@ function InitialNBody15(T)
         parse(BigFloat, "-0.00067904196080291327"),
         parse(BigFloat, "0.00032220737349778078"),
         parse(BigFloat, "-0.00314357639364532859"),
-        parse(BigFloat, "-0.00107794975959731297"),
+        parse(BigFloat, "-0.00107794975959731297")
     ]
 
     # Ceres, Pallas, Vesta, Iris, Bamberga
@@ -98,54 +98,54 @@ function InitialNBody15(T)
     qCeres = qsun + [
         parse(BigFloat, "1.438681809676469747"),
         parse(BigFloat, "-2.204373633189407045"),
-        parse(BigFloat, "-1.326397853361325874"),
+        parse(BigFloat, "-1.326397853361325874")
     ]
     qPallas = qsun + [
         parse(BigFloat, "0.203832272462290465"),
         parse(BigFloat, "-3.209619436062307152"),
-        parse(BigFloat, "0.623843179079393351"),
+        parse(BigFloat, "0.623843179079393351")
     ]
     qVesta = qsun + [
         parse(BigFloat, "0.182371836377417107"),
         parse(BigFloat, "2.386628211277654010"),
-        parse(BigFloat, "0.924596062836265498"),
+        parse(BigFloat, "0.924596062836265498")
     ]
     qIris = qsun + [
         parse(BigFloat, "1.892475267790300286"),
         parse(BigFloat, "-0.848414748075139946"),
-        parse(BigFloat, "-0.157159319044464590"),
+        parse(BigFloat, "-0.157159319044464590")
     ]
     qBamberga = qsun + [
         parse(BigFloat, "1.398759064223541682"),
         parse(BigFloat, "-1.287476729008325105"),
-        parse(BigFloat, "-0.669098428660833799"),
+        parse(BigFloat, "-0.669098428660833799")
     ]
 
     vsun = v[1:3]
     vCeres = vsun + [
         parse(BigFloat, "0.008465406136316316"),
         parse(BigFloat, "0.004684247977335608"),
-        parse(BigFloat, "0.000466157738595739"),
+        parse(BigFloat, "0.000466157738595739")
     ]
     vPallas = vsun + [
         parse(BigFloat, "0.008534313855651248"),
         parse(BigFloat, "-0.000860659210123161"),
-        parse(BigFloat, "-0.000392901992572746"),
+        parse(BigFloat, "-0.000392901992572746")
     ]
     vVesta = vsun + [
         parse(BigFloat, "-0.010174496747119257"),
         parse(BigFloat, "0.000041478190529952"),
-        parse(BigFloat, "0.001344157634155624"),
+        parse(BigFloat, "0.001344157634155624")
     ]
     vIris = vsun + [
         parse(BigFloat, "0.002786950314570632"),
         parse(BigFloat, "0.011314057384917047"),
-        parse(BigFloat, "0.004975132577079665"),
+        parse(BigFloat, "0.004975132577079665")
     ]
     vBamberga = vsun + [
         parse(BigFloat, "0.007164363244556328"),
         parse(BigFloat, "0.009219958777618218"),
-        parse(BigFloat, "0.006857861727407507"),
+        parse(BigFloat, "0.006857861727407507")
     ]
 
     q = [q; qCeres; qPallas; qVesta; qIris; qBamberga]
diff --git a/ODEProblems/InitialNBody16.jl b/ODEProblems/InitialNBody16.jl
index 8778cb5..0923b36 100644
--- a/ODEProblems/InitialNBody16.jl
+++ b/ODEProblems/InitialNBody16.jl
@@ -19,7 +19,7 @@ function InitialNBody16(T = Float64)
         parse(BigFloat, "0.310444819893871300e-13"),
         parse(BigFloat, "0.385475018780881000e-13"),
         parse(BigFloat, "0.213643444257140700e-14"),
-        parse(BigFloat, "0.138862658985619900e-14"),
+        parse(BigFloat, "0.138862658985619900e-14")
     ]
 
     N = length(Gm)
@@ -54,7 +54,7 @@ function InitialNBody16(T = Float64)
         parse(BigFloat, "-9.40015796880239402236"),
         parse(BigFloat, "-30.48331376718383944523"),
         parse(BigFloat, "-0.87240555684104999295"),
-        parse(BigFloat, "8.91157617249954997764"),
+        parse(BigFloat, "8.91157617249954997764")
     ]
 
     v = [
@@ -87,18 +87,18 @@ function InitialNBody16(T = Float64)
         parse(BigFloat, "-0.00067904196080291327"),
         parse(BigFloat, "0.00032220737349778078"),
         parse(BigFloat, "-0.00314357639364532859"),
-        parse(BigFloat, "-0.00107794975959731297"),
+        parse(BigFloat, "-0.00107794975959731297")
     ]
 
     dqMoon = [
         parse(BigFloat, "-0.00080817735147818490"),
         parse(BigFloat, "-0.0019946299854970130"),
-        parse(BigFloat, "-0.00108726268307068900"),
+        parse(BigFloat, "-0.00108726268307068900")
     ]
     dvMoon = [
         parse(BigFloat, "0.00060108481561422370"),
         parse(BigFloat, "-0.00016744546915764980"),
-        parse(BigFloat, "-0.00008556214140094871"),
+        parse(BigFloat, "-0.00008556214140094871")
     ]
 
     qq = copy(q)
@@ -119,54 +119,54 @@ function InitialNBody16(T = Float64)
     qCeres = qsun + [
         parse(BigFloat, "1.438681809676469747"),
         parse(BigFloat, "-2.204373633189407045"),
-        parse(BigFloat, "-1.326397853361325874"),
+        parse(BigFloat, "-1.326397853361325874")
     ]
     qPallas = qsun + [
         parse(BigFloat, "0.203832272462290465"),
         parse(BigFloat, "-3.209619436062307152"),
-        parse(BigFloat, "0.623843179079393351"),
+        parse(BigFloat, "0.623843179079393351")
     ]
     qVesta = qsun + [
         parse(BigFloat, "0.182371836377417107"),
         parse(BigFloat, "2.386628211277654010"),
-        parse(BigFloat, "0.924596062836265498"),
+        parse(BigFloat, "0.924596062836265498")
     ]
     qIris = qsun + [
         parse(BigFloat, "1.892475267790300286"),
         parse(BigFloat, "-0.848414748075139946"),
-        parse(BigFloat, "-0.157159319044464590"),
+        parse(BigFloat, "-0.157159319044464590")
     ]
     qBamberga = qsun + [
         parse(BigFloat, "1.398759064223541682"),
         parse(BigFloat, "-1.287476729008325105"),
-        parse(BigFloat, "-0.669098428660833799"),
+        parse(BigFloat, "-0.669098428660833799")
     ]
 
     vsun = vv[1:3]
     vCeres = vsun + [
         parse(BigFloat, "0.008465406136316316"),
         parse(BigFloat, "0.004684247977335608"),
-        parse(BigFloat, "0.000466157738595739"),
+        parse(BigFloat, "0.000466157738595739")
     ]
     vPallas = vsun + [
         parse(BigFloat, "0.008534313855651248"),
         parse(BigFloat, "-0.000860659210123161"),
-        parse(BigFloat, "-0.000392901992572746"),
+        parse(BigFloat, "-0.000392901992572746")
     ]
     vVesta = vsun + [
         parse(BigFloat, "-0.010174496747119257"),
         parse(BigFloat, "0.000041478190529952"),
-        parse(BigFloat, "0.001344157634155624"),
+        parse(BigFloat, "0.001344157634155624")
     ]
     vIris = vsun + [
         parse(BigFloat, "0.002786950314570632"),
         parse(BigFloat, "0.011314057384917047"),
-        parse(BigFloat, "0.004975132577079665"),
+        parse(BigFloat, "0.004975132577079665")
     ]
     vBamberga = vsun + [
         parse(BigFloat, "0.007164363244556328"),
         parse(BigFloat, "0.009219958777618218"),
-        parse(BigFloat, "0.006857861727407507"),
+        parse(BigFloat, "0.006857861727407507")
     ]
 
     qq = [qq; qMoon; qCeres; qPallas; qVesta; qIris; qBamberga]
diff --git a/ODEProblems/InitialNBody9.jl b/ODEProblems/InitialNBody9.jl
index d79efea..f648f99 100644
--- a/ODEProblems/InitialNBody9.jl
+++ b/ODEProblems/InitialNBody9.jl
@@ -15,7 +15,7 @@ function InitialNBody9(T)
         parse(BigFloat, "0.845970607324503000e-7"),
         parse(BigFloat, "0.129202482578296000e-7"),
         parse(BigFloat, "0.152435734788511000e-7"),
-        parse(BigFloat, "0.217844105197418000e-11"),
+        parse(BigFloat, "0.217844105197418000e-11")
     ]
 
     N = length(Gm)
@@ -50,7 +50,7 @@ function InitialNBody9(T)
         parse(BigFloat, "-9.40015796880239402236"),
         parse(BigFloat, "-30.48331376718383944523"),
         parse(BigFloat, "-0.87240555684104999295"),
-        parse(BigFloat, "8.91157617249954997764"),
+        parse(BigFloat, "8.91157617249954997764")
     ]
 
     v = [
@@ -83,7 +83,7 @@ function InitialNBody9(T)
         parse(BigFloat, "-0.00067904196080291327"),
         parse(BigFloat, "0.00032220737349778078"),
         parse(BigFloat, "-0.00314357639364532859"),
-        parse(BigFloat, "-0.00107794975959731297"),
+        parse(BigFloat, "-0.00107794975959731297")
     ]
 
     q0 = reshape(q, 3, :)
diff --git a/ODEProblems/InitialNLS.jl b/ODEProblems/InitialNLS.jl
index 2a76b27..032fb61 100644
--- a/ODEProblems/InitialNLS.jl
+++ b/ODEProblems/InitialNLS.jl
@@ -15,7 +15,7 @@ function InitialNLS(T = Float64)
         parse(BigFloat, "0.01"),
         parse(BigFloat, "0.01"),
         parse(BigFloat, "0.01"),
-        parse(BigFloat, "0.01"),
+        parse(BigFloat, "0.01")
     ]
 
     p = [
@@ -23,7 +23,7 @@ function InitialNLS(T = Float64)
         parse(BigFloat, "0.0"),
         parse(BigFloat, "0.0"),
         parse(BigFloat, "0.0"),
-        parse(BigFloat, "0.0"),
+        parse(BigFloat, "0.0")
     ]
 
     u0 = convert.(T, [q; p])
diff --git a/ODEProblems/InitialPleiades.jl b/ODEProblems/InitialPleiades.jl
index 3c48557..5ea1622 100644
--- a/ODEProblems/InitialPleiades.jl
+++ b/ODEProblems/InitialPleiades.jl
@@ -14,7 +14,7 @@ function InitialPleiades(T)
         parse(BigFloat, "0.0"),
         parse(BigFloat, "0.0"),
         parse(BigFloat, "-4.0"),
-        parse(BigFloat, "4.0"),
+        parse(BigFloat, "4.0")
     ]
 
     v = [
@@ -31,7 +31,7 @@ function InitialPleiades(T)
         parse(BigFloat, "-1.25"),
         parse(BigFloat, "1.0"),
         parse(BigFloat, "0.0"),
-        parse(BigFloat, "0.0"),
+        parse(BigFloat, "0.0")
     ]
 
     u0 = Array{T}(undef, 1)
diff --git a/src/IRKCoefficients.jl b/src/IRKCoefficients.jl
index 88d78a6..30b16aa 100644
--- a/src/IRKCoefficients.jl
+++ b/src/IRKCoefficients.jl
@@ -164,7 +164,7 @@ function HCoefficients!(mu, hc, hb, nu, h, hprev, T::Type{<:CompiledFloats})
             +1.8134189168918099148257522463859781e-01,
             +1.5685332293894364366898110099330067e-01,
             +1.1119051722668723527217799721312045e-01,
-            +5.0614268145188129576265677154981094e-02,
+            +5.0614268145188129576265677154981094e-02
         ])
 
     hc .= convert.(T,
@@ -176,7 +176,7 @@ function HCoefficients!(mu, hc, hb, nu, h, hprev, T::Type{<:CompiledFloats})
             +5.9171732124782490246973807118009203e-01,
             +7.6276620495816449290886952459462321e-01,
             +8.9833323870681336979577696823791522e-01,
-            +9.8014492824876811584178043428473653e-01,
+            +9.8014492824876811584178043428473653e-01
         ])
 
     s = length(hb)
@@ -217,7 +217,7 @@ function HCoefficients!(mu, hc, hb, nu, h, hprev, T)
         parse(T, "+1.8134189168918099148257522463859781e-01"),
         parse(T, "+1.5685332293894364366898110099330067e-01"),
         parse(T, "+1.1119051722668723527217799721312045e-01"),
-        parse(T, "+5.0614268145188129576265677154981094e-02"),
+        parse(T, "+5.0614268145188129576265677154981094e-02")
     ]
 
     hc .= [
@@ -228,7 +228,7 @@ function HCoefficients!(mu, hc, hb, nu, h, hprev, T)
         parse(T, "+5.9171732124782490246973807118009203e-01"),
         parse(T, "+7.6276620495816449290886952459462321e-01"),
         parse(T, "+8.9833323870681336979577696823791522e-01"),
-        parse(T, "+9.8014492824876811584178043428473653e-01"),
+        parse(T, "+9.8014492824876811584178043428473653e-01")
     ]
 
     s = length(hb)
@@ -267,7 +267,7 @@ function EstimateCoeffs!(alpha, T::Type{<:CompiledFloats})
             +5.9171732124782490246973807118009203e-01,
             +7.6276620495816449290886952459462321e-01,
             +8.9833323870681336979577696823791522e-01,
-            +9.8014492824876811584178043428473653e-01,
+            +9.8014492824876811584178043428473653e-01
         ])
 
     s = length(c)
@@ -291,7 +291,7 @@ function EstimateCoeffs!(alpha, T)
         parse(T, "+5.9171732124782490246973807118009203e-01"),
         parse(T, "+7.6276620495816449290886952459462321e-01"),
         parse(T, "+8.9833323870681336979577696823791522e-01"),
-        parse(T, "+9.8014492824876811584178043428473653e-01"),
+        parse(T, "+9.8014492824876811584178043428473653e-01")
     ]
 
     s = length(c)
diff --git a/src/IRKGL16AuxFunctions.jl b/src/IRKGL16AuxFunctions.jl
index a11592b..62832b6 100644
--- a/src/IRKGL16AuxFunctions.jl
+++ b/src/IRKGL16AuxFunctions.jl
@@ -26,6 +26,10 @@ function ErrorEst(U, F, dt, beta, abstol, reltol)
         end
     end
 
+    if est == 0
+        est = eps(realuiType)
+    end
+
     return (est / D)
 end
 
diff --git a/src/IRKGL16Solver.jl b/src/IRKGL16Solver.jl
index 005fc5a..5894591 100644
--- a/src/IRKGL16Solver.jl
+++ b/src/IRKGL16Solver.jl
@@ -58,7 +58,7 @@ abstract type IRKAlgorithm{
     threading,
     mixed_precision,
     low_prec_type,
-    nrmbits,
+    nrmbits
 } <: OrdinaryDiffEqAlgorithm end
 struct IRKGL16{
     mstep,
@@ -68,7 +68,7 @@ struct IRKGL16{
     threading,
     mixed_precision,
     low_prec_type,
-    nrmbits,
+    nrmbits
 } <: IRKAlgorithm{
     mstep,
     maxtrials,
@@ -77,7 +77,7 @@ struct IRKGL16{
     threading,
     mixed_precision,
     low_prec_type,
-    nrmbits,
+    nrmbits
 } end
 function IRKGL16(;
         mstep = 1,
@@ -96,14 +96,15 @@ function IRKGL16(;
         threading,
         mixed_precision,
         low_prec_type,
-        nrmbits,
+        nrmbits
     }()
 end
 
-function DiffEqBase.__solve(prob::DiffEqBase.AbstractODEProblem{
+function DiffEqBase.__solve(
+        prob::DiffEqBase.AbstractODEProblem{
             uType,
             tspanType,
-            isinplace,
+            isinplace
         },
         alg::IRKGL16{
             mstep,
@@ -113,7 +114,7 @@ function DiffEqBase.__solve(prob::DiffEqBase.AbstractODEProblem{
             threading,
             mixed_precision,
             low_prec_type,
-            nrmbits,
+            nrmbits
         },
         args...;
         dt = 0.0,
@@ -133,8 +134,8 @@ function DiffEqBase.__solve(prob::DiffEqBase.AbstractODEProblem{
         threading,
         mixed_precision,
         low_prec_type,
-        nrmbits,
-    }
+        nrmbits
+}
     s = 8
     stats = DiffEqBase.Stats(0)
 
diff --git a/test/runtests.jl b/test/runtests.jl
index fb28d5c..f399c7f 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -45,7 +45,6 @@ sim = test_convergence(dts, prob_ode_bigfloat2Dlinear, IRKGL16())
 
 # Backward integrations tests
 
-
 #Define the problem
 const g = 9.81
 L = 1.0