using DynamicalSystems
using Plots

function σ(x,b,a)
	return 1. / (1+exp(a*(b-x)))
end

function f(u,p,t)
	# parameters
	ψ = p[1]; wS = p[2]; wy = p[3]; τx = p[4]; τb = p[5]; τS = p[6]; a = p[7]
	# aux vars
	y = σ(u[1],u[2],a)
	# equations
	du1 = (-u[1] + wS*ψ*u[3] + wy*(1-ψ)*y)/τx
	du2 = (y - 1/2)/τb
	du3 = (y - u[3])/τS
	return SVector{3}(du1,du2,du3)
end

begin
	Ψ = 0.5; wS=2; wy=2; τx = 0.2; τb = 2.; τS = 1; a = 4.
	ds = ContinuousDynamicalSystem(f, rand(3), [Ψ, wS, wy, τx, τb, τS, a])
end

traj = trajectory(ds, 1)

begin
	plot(traj[:,1],traj[:,2])
end