Plotting 3D graphic for a differential system with Manipulate

I try to adapt the following code with three differential equations to a system with four differential equations. My aim is to see a bifurcation in this differential system by looking at three variables of the system. (Mathematically, I proved the possibility of bifurcations for some exogenous parameter set)

I have the following code :

Manipulate[
With[{sol =
First[NDSolve[{{λ'[
t] == (((ρ + θ) -
a Subscript[α, 1] k^(
Subscript[α, 1] - 1)) λ[t] - μ[
t] γ a Subscript[α, 1] k^(
Subscript[α, 1] -
1) - θ (((1 - E^(
n (-ζ + θ + ρ))) (k^(-1 - α) \
α ψ + ((-1 + E^((
n (-ζ + θ + ρ))/(-1 + β + \
σ))) (a k^Subscript[α, 1])^((
1 - σ)/β) (-1 + β + σ) \
(((-1 + E^((
n β (-ζ + θ + ρ))/(-1 + \
β + σ))) (-1 + β + σ))/(-ζ + θ \
+ ρ))^(-1 + σ) Subscript[α, 1])/(
k β (-ζ + θ + ρ))))/(ζ \
- θ - ρ))), μ'[
t] == ((ρ + θ) μ[t] - μ [t] (1 - 2 s)),
k'[t] == (
a k^Subscript[α,
1] - ((λ [
t] β )/((((β + σ - 1) (1 -
Exp[-(((ζ - (ρ + θ)) n)/(β + \
σ -
1))]))/(ζ - (ρ + \
θ)))/((((β + σ - 1) (1 -
Exp[-(((ζ - (ρ + θ)) β n)/(\
β + σ - 1))]))/(ζ - (ρ + θ)))^(
1 - σ))) )^(β/(1 - β - σ))),
s'[t] == (s (1 - s) - γ a k^Subscript[α, 1])},
k[0] == 8424.78495705205 (1 - 10^-15),
s[0] == 0.45999999999999996 (1 - 10^-15), μ[
0] == -4787.8954001202255 (1 - 10^-15), λ[0] ==
96.96238348579865 (1 - 10^-15)}, {λ, μ, k, s}, {t,
0, tf}, MaxSteps -> ∞]]},
Column[{
ParametricPlot3D[
Evaluate[{k[t], λ[t], s[t]} /. sol], {t, 0, tf},
BoxRatios -> 1, PlotPoints -> 1000, PlotRange -> All,
ImageSize -> {400, 400}, SphericalRegion -> True, Ticks -> False],
Plot[Evaluate[{k[t], λ[t], s[t]} /. sol], {t, 0, tf},
PlotStyle -> Automatic, ImageSize -> 400, AspectRatio -> 1/6]}]],
{{tf, 130, Style["t", Italic]}, 1, 130, ImageSize -> Tiny},
{{θ, 0.95}, -.1, 1, ImageSize -> Tiny},
{{β, 0.65}, 0.18, 3, ImageSize -> Tiny},
{{ρ, 1}, 0.01, 1, ImageSize -> Tiny},
{{Subscript[α, 1], 1}, -.1, 1, ImageSize -> Tiny},
{{ψ, 100}, -.1, 100, ImageSize -> Tiny},
{{γ, 0.05}, -.0 .01, 0.05, ImageSize -> Tiny},
{{ζ, 0.05}, -.0 .01, 0.05, ImageSize -> Tiny},
{{σ, 0.05}, -.0 .1, 2, ImageSize -> Tiny},
{{α, 0.05}, -.0 .1, 2, ImageSize -> Tiny},
{{a, 0.05}, -.0 .1, 4, ImageSize -> Tiny},
{{n, 90}, -.0 .1, 90, ImageSize -> Tiny},
ControlPlacement -> Left, SynchronousUpdating -> False]


I know that equations look horrible, sorr for that. When I remove NDSolve part from Manipulate and attribute all parameters by numerical values below (there are 11 constant parameters in total + time $t$), everythings work but I don't know why the code does not give any output with Manipulate. What am I missing ?

paramFinal = {ρ -> 0.05, θ -> 0.03, n -> 80, Subscript[α, 1] -> 0.3, α -> 0.65, β -> 1.05, a -> 1.65, γ -> 0.01, σ -> 0.5, ψ -> 1, ζ -> 0.014870264736785313};

• A hint: extract the NDSolve part from your Manipulate and fix ts to 130. You will see that your parameter and function substitutions contain errors. The final equations you use in NDSolve still contain unsubstituted parameters and you don't get a solution. When you fixed this and get a solution, you can take care of your Manipulate but not sooner. – halirutan Feb 25 '17 at 3:56
• How can you vary d as a parameter, whereas d is defined as an expression? – zhk Feb 25 '17 at 6:04
• @halirutan Thanks for the hint. When I extract the NDSolve from Manipulate as you have mentionned and attribute numerical values to parameters, everything works fine but with Manipulate, there are no error messages but any output neither. – optimal control Feb 25 '17 at 14:19
• @optimalcontrol I don't see how that's possible. You have multiple instances in your code of k and s appearing without their argument, for which NDSolve complains loudly. Are you sure you are using the exact same code in your notebook? – MarcoB Feb 26 '17 at 0:36