function renormHW1_20060919
%This is the matlab code for HW1, by Ed Greco

%HW Problem 26.1-----
%--------------------

%B) discretize lambda and re-iterate using the last x_n as x_0
clear lambda x
lambda = linspace(.01,2,1e3);
x(1) = .5;
num_iter = 1e4;
figure
for j = 1:1e3
  for i = 1:num_iter
    x(i+1) = lambda(j) - x(i)^2;
  end
  plot(lambda(j)*ones(200,1)',x(num_iter+1-199:end),'.')
  hold on
  x(1) = x(end);
end
hold off
xlabel('lambda')
ylabel('asymptotic fixed points')


%C) the critical point is x = 0. for x_0 = 0,
%   x_1 = lambda.  If x_1 is greater than the negative 
%   of the negative fixed point, it will iterate to -infinity.
%   f(x) = x => x'^2+x'-lambda = 0 then
%   -x_' = 1/2 + sqrt(1/4+lambda)
%   and 1/2 + sqrt(1/4+lambda) = lambda => lambda = 2.  

%D) try lambda = 3/4
clear lambda x
lambda = 3/4;
x(1) = .25;
num_iter = 1e4; %the number of times to iterate
for i = 1:num_iter
  x(i+1) = lambda - x(i)^2;
end
figure
plot([1:num_iter+1],x,'.')
xlabel('n')
ylabel('x_n')
% The convergance is much slower...

%E) f(x) = x => x'^2+x'-lambda = 0 then
%   x' = -1/2 +- sqrt(1+4*lambda)
%   for stability abs( df/dx ) < 1
%     abs(-2*x') < 1 then
%     -1 < 1 -+ sqrt(1+4*lambda) < 1
%      4 > 1+4*lambda > 0
%      3/4 > lambda > -1/4
%   then for lambda = 3/4 the fixed point is unstable

%F) try differnt lambda
clear lambda x
lambda = [1,1.31070274134,1.38154748443,1.3979453597];
x(1) = .25;
num_iter = 1e4; %the number of times to iterate
figure
for j = 1:4
  for i = 1:num_iter
    x(i+1) = lambda(j) - x(i)^2;
  end
  subplot(4,1,j)
  plot([1:num_iter+1],x,'.')
  xlabel('n')
  ylabel('x_n')
  x(1) = x(end);
end

%G) to do by newtons method...
%H) to do...

%HW Problem 26.2-----
%--------------------

%A) dx_1 = x_2
%   dx_2 = A*cos(w*t) - k*x_2 + x_1 - 4*x_1^3

%B) dx_1 = x_2
%   dx_2 = x_3 - k*x_2 + x_1 - 4*x_1^3
%   dx_3 = x_4
%   dx_4 = -x_3/w^2

%C)  
figure
k = 0.154;
[Y] = solveode(.1,k,[0 250*(2*pi/1.2199778)],[0.5,0]);
plot(Y(:,1),Y(:,2),'.')
xlabel('x')
ylabel('dx/dt')
clear Y

%D) See attached plots
[Y] = solveode(.1,k,[0 250*2*pi/1.2199778],[0.5 0]);
x_start = [Y(end,1),Y(end,2)];
clear Y
figure
for j = 1:9
  if j == 1
    [Y] = solveode(.1,k,[0 (2*pi/1.2199778)],x_start);
  else
    [Y] = solveode(.1,k,[2^(j-2)*(2*pi/1.2199778) 2^(j-1)*(2*pi/1.2199778)],x_start);
  end
  pltd(j,1) = Y(end,1);
  pltd(j,2) = Y(end,2);
  x_start = [Y(end,1),Y(end,2)];
  clear Y
end
plot(pltd(:,1),pltd(:,2),'.')

%E) Lets try plotting the return map for the previous Poincare section
figure
subplot(2,2,1)
plot(pltd(1:(end-1),1),pltd(2:end,1),'.')
xlabel('x_n')
ylabel('x_{n+1}')
subplot(2,2,2)
plot(pltd(1:(end-1),1),pltd(2:end,2),'.')
xlabel('x_n')
ylabel('dx_{n+1}')
subplot(2,2,3)
plot(pltd(1:(end-1),2),pltd(2:end,1),'.')
xlabel('dx_n')
ylabel('x_{n+1}')
subplot(2,2,4)
plot(pltd(1:(end-1),2),pltd(2:end,2),'.')
xlabel('dx_n')
ylabel('dx_{n+1}')
clear pltd

%F) period doubling is so cool, too bad it doesnt work for me in this case.
clear A
A = linspace(.05,.2,1e2);
[Y] = solveode(A(1),k,[0 250*2*pi/1.2199778],[0.5 0]);
x_start = [Y(end,1),Y(end,2)];
figure
for i = 1:1e2
  for j = 1:10
    [Y] = solveode(A(i),k,[(j-1)*(2*pi/1.2199778) j*(2*pi/1.2199778)],x_start);
    x_start = [Y(end,1),Y(end,2)];
    plot(A(i)*ones(200,1)',Y(end-199:end,1),'.')
    hold on
    clear Y
  end  
end
clear A

%G) See attached plots, almost works but some points are not bifurcating?
A = [0.1 0.11 0.114 0.11437];
[Y] = solveode(A(1),k,[0 500*2*pi/1.2199778],[0.5 0]);
x_start = [Y(end,1),Y(end,2)];
clear Y 
figure
for i = 1:4
  subplot(2,2,i)
  for j = 1:10
    [Y] = solveode(A(i),k,[(j-1)*(2*pi/1.2199778) j*(2*pi/1.2199778)],x_start);
    pltd(j,1) = Y(end,1);
    pltd(j,2) = Y(end,2);
    x_start = [Y(end,1),Y(end,2)];
    clear Y
  end
  plot(pltd(:,1),pltd(:,2),'.')
  clear pltd
end

end

function [Y] = solveode(A,k,t_domain,x_start)
%Solve the ode with the given parameters for A and k, this must be done
%nested to avoid use of global assignments
  options = odeset('RelTol',1e-8,'AbsTol',1e-12);
  [T,Y] = ode45(@nonlinear_oscillator,t_domain,x_start,options);
  function dx = nonlinear_oscillator(t,x)
    %The ode from problem 27.2
    dx = zeros(2,1);    % a column vector
    dx(1) = x(2);
    dx(2) = A*cos(1.2199778*t) - k*x(2) + x(1) - 4*x(1)^3;
  end
end