#=
Test function 

Properties 
----------

- There are 9 steady states:
    - Stable
        - (-2.39,2.35)    (2.3,-2.34) (-2.18,-2.12)   (2.06,2.12) 
    - Unstable
        - (-2.27,-0.23)   (2.17,0.22) (-0.12,-2.23)   (0.33,2.22) (0.1,0.01)

Reference 
---------
1. Source code from https://static-content.springer.com/esm/art%3A10.1038%2Fs41598-017-15889-2/MediaObjects/41598_2017_15889_MOESM2_ESM.rar
=#

module SampleDE

"""
Reference
---------
- `ODE.m`
"""
function ODE!(du,u,p,t)
    du[1]= -1 + 9 *u[1] - 2*u[1]^3 + 9 *u[2] - 2*u[2]^3;
    du[2]=  1 - 11*u[1] + 2*u[1]^3 + 11*u[2] - 2*u[2]^3;
end

function ODE(u,p,t)
    du = deepcopy(u)
    return ODE!(du,u,p,t)
end

u0 = [1.0, 2.0]
p = [2.0, 3.0]



"""
Reference
---------
- `exact_potential.m`
"""
function exact_potential(points, RefPointNum)
    pot = zeros(length(points), 1)
    for i in 1:length(points)
        x = points[i][1]; y = points[i][2];
        pot[i] = -5*(x^2+y^2)+(x^4+y^4)/2+x*y+x;
    end
    pot = pot .- pot[RefPointNum]
    return pot
end


"""
Reference
---------
- `Readme.pdf`
"""
ODE_features = Dict(
    "StableStates" => [ (-2.39,2.35),(2.3,-2.34), (-2.18,-2.12),   (2.06,2.12)  ],
    "UnstableStates" => [(-2.27,-0.23) ,  (2.17,0.22), (-0.12,-2.23),   (0.33,2.22), (0.1,0.01)],
    "SteadyStates" => [(-2.39,2.35),(2.3,-2.34), (-2.18,-2.12),   (2.06,2.12), (-2.27,-0.23) ,  (2.17,0.22), (-0.12,-2.23),   (0.33,2.22) ,(0.1,0.01)]
)



function gradient(dt,x) 

    x = reshape(x, 2, 3)

    F=[ [(dt*(1/dt - 6*x[1, 1]^2 + 9)*(9*x[1, 1] + 9*x[2, 1] + (x[1, 1] - x[1, 2])/dt - 2*x[1, 1]^3 - 2*x[2, 1]^3 - 1))/2 + (dt*(6*x[1, 1]^2 - 11)*(11*x[2, 1] - 11*x[1, 1] + (x[2, 1] - x[2, 2])/dt + 2*x[1, 1]^3 - 2*x[2, 1]^3 + 1))/2,x[1, 1]^3 - (9*x[2, 1])/2 - (x[1, 1] - x[1, 2])/(2*dt) - (9*x[1, 1])/2 + x[2, 1]^3 + (dt*(1/dt - 6*x[1, 2]^2 + 9)*(9*x[1, 2] + 9*x[2, 2] + (x[1, 2] - x[1, 3])/dt - 2*x[1, 2]^3 - 2*x[2, 2]^3 - 1))/2 + (dt*(6*x[1, 2]^2 - 11)*(11*x[2, 2] - 11*x[1, 2] + (x[2, 2] - x[2, 3])/dt + 2*x[1, 2]^3 - 2*x[2, 2]^3 + 1))/2 + 1/2,x[1, 2]^3 - (9*x[2, 2])/2 - (x[1, 2] - x[1, 3])/(2*dt) - (9*x[1, 2])/2 + x[2, 2]^3 + 1/2 ],[
    (dt*(1/dt - 6*x[2, 1]^2 + 11)*(11*x[2, 1] - 11*x[1, 1] + (x[2, 1] - x[2, 2])/dt + 2*x[1, 1]^3 - 2*x[2, 1]^3 + 1))/2 - (dt*(6*x[2, 1]^2 - 9)*(9*x[1, 1] + 9*x[2, 1] + (x[1, 1] - x[1, 2])/dt - 2*x[1, 1]^3 - 2*x[2, 1]^3 - 1))/2,(11*x[1, 1])/2 - (11*x[2, 1])/2 - (x[2, 1] - x[2, 2])/(2*dt) - x[1, 1]^3 + x[2, 1]^3 + (dt*(1/dt - 6*x[2, 2]^2 + 11)*(11*x[2, 2] - 11*x[1, 2] + (x[2, 2] - x[2, 3])/dt + 2*x[1, 2]^3 - 2*x[2, 2]^3 + 1))/2 - (dt*(6*x[2, 2]^2 - 9)*(9*x[1, 2] + 9*x[2, 2] + (x[1, 2] - x[1, 3])/dt - 2*x[1, 2]^3 - 2*x[2, 2]^3 - 1))/2 - 1/2,(11*x[1, 2])/2 - (11*x[2, 2])/2 - (x[2, 2] - x[2, 3])/(2*dt) - x[1, 2]^3 + x[2, 2]^3 - 1/2]]
    return hcat(F...)'
end


end