clearconsole()
using Plots
using LinearAlgebra
using DifferentialEquations
using Measures

##PlotVectorField
function PlotVectorField(vectorx1,vectorx2,epsilon,factor)
    # uncomment next four lines if you want to plot the tangent vectors for a
    # region of the state-plane
    #Nx1=length(vectorx1)
    #Nx2=length(vectorx2)
    #xx1=repeat(vectorx1,Nx2)
    #xx2=repeat(vectorx2,inner=Nx1)
    xx1=vectorx1
    xx2=vectorx2
    u=xx2
    v=-xx1-epsilon*((xx1.^2).-1).*xx2
    uplot=u./sqrt.(u.^2+v.^2)
    vplot=v./sqrt.(u.^2+v.^2)
    fig=quiver(xx1,xx2,quiver=(uplot/factor,vplot/factor),
    xlims=(minimum(vectorx1),maximum(vectorx1)),
    ylim=(minimum(vectorx2),maximum(vectorx2)),aspect_ratio=1,color=:black,linecolor=:black)
    return fig
end

# Integrate the van der pol equation
function vanderPol!(dx,x,epsilon,t)
    dx[1] = x[2]
    dx[2] = -x[1]-epsilon*(x[1]^2-1)*x[2]
end

## Parameters
epsilon=0.3
x1min=-4.0
x1max=4.0
x2min=-4.0
x2max=4.0
dx1=.6
dx2=.6
factor=1.5
Tmax=1000
u0=0.0
u0dot=0.1

vectorx1=collect(x1min:dx1:x1max) # Employed if you want to plot the vector field
                                  # for a region of the phase-space
vectorx2=collect(x2min:dx2:x2max) # Employed if you want to plot the vector field
                                  # for a region of the phase-space
x01=[u0;u0dot] ## First initial condition to be simulated
x02=[2.5;-2.5]  ## Second initial condition to be simulated

# Simulations
tspan = (0.0,Tmax)
prob = ODEProblem(vanderPol!,x01,tspan,epsilon)
sol1 = solve(prob)
prob2 = ODEProblem(vanderPol!,x02,tspan,epsilon)
sol2 = solve(prob2)

tArrow=collect(0.0:Tmax/20:Tmax) #Instants in which the tangent vector are
                                 #calculated

# Manipulating the solutions
sol1ArrayArrow=Transpose(Array(sol1(tArrow)))
sol1uArrow=sol1ArrayArrow[:,1]
sol1upArrow=sol1ArrayArrow[:,2]

sol2ArrayArrow=Transpose(Array(sol2(tArrow)))
sol2uArrow=sol2ArrayArrow[:,1]
sol2upArrow=sol2ArrayArrow[:,2]

# Plotting
plotu=plot(sol1,vars=(1),linewidth=1.5,xtickfontsize=14,ytickfontsize=14,
xlims=(0.0,Tmax),ylims=(-2.5,2.5),linecolor=:orangered,legend=false,
right_margin=5mm,left_margin=-2mm,framestyle=:box)
xlabel!("t",xguidefontsize=16)
ylabel!("x₁",yguidefontsize=16)
display(plotu)
NameFig="vanderPolx1.pdf"
savefig(NameFig)


fig=PlotVectorField(sol1uArrow,sol1upArrow,epsilon,factor)
plot!(fig,sol1,vars=(1,2),linewidth=1.5,xtickfontsize=14,ytickfontsize=14,
xlims=(-2.5,2.5),ylims=(-2.5,2.5),linecolor=:orangered,legend=false,
right_margin=5mm,left_margin=-2mm,framestyle=:box)
xlabel!("x₁",xguidefontsize=16)
ylabel!("x₂",yguidefontsize=16)
display(fig)
NameFig="vanderPolx1x2vec.pdf"
savefig(NameFig)


plotuup=plot(sol1,vars=(1,2),linewidth=1.5,xtickfontsize=14,ytickfontsize=14,
xlims=(-3.5,3.5),ylims=(-3.5,3.5),linecolor=:orangered,legend=false,
right_margin=5mm,left_margin=-2mm,framestyle=:box)
plot!(sol2,vars=(1,2),linewidth=1.5,xtickfontsize=14,ytickfontsize=14,
xlims=(-3.5,3.5),ylims=(-3.5,3.5),linecolor=:blue2,legend=false)
xlabel!("x₁",xguidefontsize=16)
ylabel!("x₂",yguidefontsize=16)
display(plotuup)
NameFig="vanderPolx1x22ICs.pdf"
savefig(NameFig)