make_rp;

Nlevs=rp.Nn;
n=rp.Nn;
rp.a1=-0.1950; rp.a2=0.95; rp.lam1=1.1; rp.lam2=0.9;rp.Jmin = 0.0965;
[v, w]=steepest_descent(rp);
Jn1  = rp.lam1.*(1-rp.mu.*rp.lam1).^(2*n).*v(1,1:Nlevs).^2 + rp.lam2.*(1-rp.mu.* rp.lam2).^(2*n).*v(2,1:Nlevs).^2 + rp.Jmin;

rp.a1=-0.9750; rp.a2=0.95; rp.lam1=1.5;   rp.lam2=0.5; rp.Jmin = 0.0731;
[v, w]=steepest_descent(rp);
Jn2  = rp.lam1.*(1-rp.mu.*rp.lam1).^(2*(1:n)).*v(1,1:Nlevs).^2 + rp.lam2.*(1-rp.mu*rp.lam2).^(2*(1:n)).*v(2,1:Nlevs).^2+ rp.Jmin;;

rp.a1=-1.5955; rp.a2=0.95; rp.lam1=1.818; rp.lam2=0.182; rp.Jmin = 0.0322;
[v, w]=steepest_descent(rp);
Jn3  = rp.lam1.*(1-rp.mu*rp.lam1).^(2*(1:n)).*v(1,1:Nlevs).^2 + rp.lam2.*(1-rp.mu*rp.lam2).^(2*(1:n)).*v(2,1:Nlevs).^2+ rp.Jmin;;

rp.a1=-1.9114; rp.a2=0.95; rp.lam1=1.957; rp.lam2=0.0198;rp.Jmin = 0.0038;
[v, w]=steepest_descent(rp);
Jn4  = rp.lam1.*(1-rp.mu*rp.lam1).^(2*(1:n)).*v(1,1:Nlevs).^2 + rp.lam2.*(1-rp.mu*rp.lam2).^(2*(1:n)).*v(2,1:Nlevs).^2+ rp.Jmin;;

plot(Jn1,'b-')
hold on
plot(Jn2,'r:')
plot(Jn3,'k--')
plot(Jn4,'g-.')
hold off
legend('\chi = 1.22', '\chi = 3', '\chi = 10', '\chi = 100')
xlabel('Time, n')
ylabel('J(n)')

print -dpsc 4_10