"Ramsey growth model";
"Las ecuaciones diferenciales del modelo son:";
"k'[t]==f[k[t]]-(δ+n)*k[t]-c[t]";
"c'[t]==(1/σ)*(f'[k]-(δ+β))*c[t]";

"Las vamos a calcular para los siguientes valores paramétricos y para la siguiente función de producción";
β=0.03;σ=1/2; α=0.25; n=0.01; δ=0.05;
f[k_]=(k^α);

"Isoclina de k y su representación";
Dkrespt[c_,k_]=f[k]-(δ+n)*k-c;
Isok=Solve[Dkrespt[c,k]==0,c];
Isokk[k_]=Isok[[1,1,2]];
Graphkpunto=Plot[Isokk[k],{k,0,45},PlotRange->{{0,45},{0,2}},PlotStyle->Thick]

"Isoclina de c y su representación. Tener en cuenta que es vertical";
kfi=FindRoot[(1/σ)*(f'[k]-(δ+β))==0,{k,1}];
kfijo=kfi[[1,2]]
Graphcpunto=ParametricPlot[{kfijo,c},{c,0,2},PlotRange->{{0,45},{0,2}},AspectRatio->1/GoldenRatio,PlotStyle->Thick]

"Isoclinas de k y de c y su representación conjunta";
Isoclinas=Show[{Graphkpunto,Graphcpunto}]

"Campo vectorial y las dos isoclinas";
vectores=VectorPlot[{f[k]-(δ+n)*k-c,((1/σ)*(f'[k]-(δ+β)))*c},{k,1,45},{c,0,2},AspectRatio->1/GoldenRatio,PlotRange->{{0,45},{0,2}}];
Show[{vectores,Graphkpunto,Graphcpunto},Axes->True,AxesOrigin->{0,0}]

"Punto fijo y comprobación del caracter de punto silla del punto fijo";
"El punto fijo es";
cfijo=f[k]-(δ+n)*k/.k->kfijo;
{kfijo,cfijo}

"Matriz de la aproximación lineal en el punto fijo";
a11=D[Dkrespt[c,k],k]/.k->kfijo/.c->cfijo;
a12=D[Dkrespt[c,k],c]/.k->kfijo/.c->cfijo;
Dcrespt[c_,k_]=(1/σ)(f'[k]-δ-β)*c;
a21=D[Dcrespt[c,k],k]/.k->kfijo/.c->cfijo;
a22=D[Dcrespt[c,k],c]/.k->kfijo/.c->cfijo;
A={{a11,a12},{a21,a22}};
MatrixForm[A]

"Valores y vectores propios";
"valores propios"
ValProp=Eigenvalues[A];
{ValProp[[1]],ValProp[[2]]}
"Vectores propios"
VectProp=Eigenvectors[A];
{MatrixForm[VectProp[[1]]],MatrixForm[VectProp[[2]]]}
"luego es punto silla"

"Dinamica y campo de vectores";
Derf[k_]=f'[k];
eq1:={ k'[t]==f[k[t]]-(δ+n)*k[t]-c[t],c'[t]==((1/σ)*(Derf[k[t]]-(δ+β)))*c[t],k[0]==9.2,c[0]==2}
eq2:={ k'[t]==f[k[t]]-(δ+n)*k[t]-c[t],c'[t]==((1/σ)*(Derf[k[t]]-(δ+β)))*c[t],k[0]==0.4,c[0]==0.23}
eq3:={ k'[t]==f[k[t]]-(δ+n)*k[t]-c[t],c'[t]==((1/σ)*(Derf[k[t]]-(δ+β)))*c[t],k[0]==20,c[0]==2}
eq4:={ k'[t]==f[k[t]]-(δ+n)*k[t]-c[t],c'[t]==((1/σ)*(Derf[k[t]]-(δ+β)))*c[t],k[0]==2,c[0]==0.75}
var:={k,c}
trange:={t,0,60}
sol1=NDSolve[eq1,var,trange];
sol2=NDSolve[eq2,var,trange];
sol3=NDSolve[eq3,var,trange];
sol4=NDSolve[eq4,var,trange];
graph1=ParametricPlot[Evaluate[{k[t],c[t]}/.sol1],{t,0,25},PlotPoints->500,PlotRange->{{0,45},{0,2}},PlotStyle->Thick];
graph2=ParametricPlot[Evaluate[{k[t],c[t]}/.sol2],{t,0,50},PlotPoints->500,PlotRange->{{0,45},{0,2}},PlotStyle->Thick];
graph3=ParametricPlot[Evaluate[{k[t],c[t]}/.sol3],{t,0,50},PlotPoints->500,PlotRange->{{0,45},{0,2}},PlotStyle->Thick];
graph4=ParametricPlot[Evaluate[{k[t],c[t]}/.sol4],{t,0,9},PlotPoints->500,PlotRange->{{0,45},{0,2}},PlotStyle->Thick];
Show[{Isoclinas,vectores,graph1,graph2,graph3,graph4},Axes->True,AxesOrigin->{0,0}]