%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Example Octave/MATLAB programs to produce the graphs
%% for Math 3520, HW 5, Spring 2026
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%% Ex.4
%%%% Runge-Kutta 4th order solver in 3D:
function X = rk4b(X0,F,steps,h)
 % Input
 %  X0 = [x0,y0,z0]   Row vector, initial position
 %  F  = function handle, computes Xdot(x,y,z)
 %  steps = number of time steps to compute (e.g., 100)
 %  h = time increment per step (e.g., 1/steps )
 % Output
 %  X = (1+steps)x(3) matrix of (x,y,z) values along
 %                    the trajectory starting at X0
  X = zeros(1+steps,3);  % Solution to be computed
  X(1,:) = X0;  % Initial point
  for n=1:steps
    f1 = F(X(n,:));
    f2 = F(X(n,:)+f1*h/2);
    f3 = F(X(n,:)+f2*h/2);
    f4 = F(X(n,:)+f3*h);
    X(n+1,:) = X(n,:) +h*(f1+2*f2+2*f3+f4)/6;
  end
end
%%%% Function to find local arg-maxima and maximum values:
function [Xs,Ys] = lmax(t,z)
% Input
%  t = Row vector of independent variable samples
%  z = row vector of dependent variable samples
% Output
%  X = row vector of local maxima
%  Y = row vector local maximum values
 dz = diff(z);  % differences z(k+1)-z(k), k=1,2,...
 n = 0;         % initial local max index
 for k=1:(length(dz)-1)
   if ( dz(k)>0 && dz(k+1)<0 )  % found local max at k+1
      n = n+1;      % increment the count of local maxes
      x1=t(k); x2=t(k+1); x3=t(k+2);
      y1=z(k); y2=z(k+1); y3=z(k+2);
      xnum=x1*x1*(y3-y2)+x2*x2*(y1-y3)+x3*x3*(y2-y1);
      xden=x1*(y3-y2)+x2*(y1-y3)+x3*(y2-y1);
      Xs(n) = 0.5*xnum/xden;  % append nth arg max
      d12=x1-x2; d31=x3-x1; d23=x2-x3;
      ynum1= 2*(y1*y1*d23^4+y2*y2*d31^4+y3*y3*d12^4);
      ynum2= (y1*d23^2+y2*d31^2+y3*d12^2)^2;
      yden = 4*d12*d23*d31*(y1*d23+y2*d31+y3*d12);
      Ys(n) = (ynum1-ynum2)/yden; % append nth max val
   end  % endif
 end  % endfor
 return
end % endfunction
% Compute the Lorenz map from a particular trajectory
%%%% From Strogatz, 3rd edition, p.353
X0=[0,1,0]; % Lorenz's initial value
sigm=10; % Prandtl number, from Lorenz
rayl=28; % Rayleigh number, from Lorenz
b=8/3;   % unnamed aspect ratio, from Lorenz
Xdot = @(V) [sigm*(V(2)-V(1)), (rayl-V(3))*V(1)-V(2), V(1)*V(2)-b*V(3)];
X = rk4b(X0, Xdot, 10000,h);  % Complete trajectory to t=100.
%% Plot the Z coordinate:
t = (0:10000)*h;     % time values
z = X(:,3);          % z-coordinate values
figure; plot(t,z);   % display the trajectory's z-coordinate
title("Z-coordinate of Lorenz System Trajectory, r=28")
xlabel("T"); ylabel("Z(T)");
%% Plot the Lorenz map:
[Xs,Ys] = lmax(t,z); % find the local maximum values in z()
Zn = Ys(1:length(Ys)-1);
Znp1 = Ys(2:length(Ys));
figure;  plot(Zn,Znp1,"o");  % point plot of z_n vs. z_{n+1}
hold on;
zz=floor(min(Ys)):ceil(max(Ys)); plot(zz,zz);  % plot Y=X line
title("Lorenz map, r=28"); xlabel("Z_n"); ylabel("Z_{n+1}");

%%% Ex.5
%%%% From Strogatz, 3rd edition, p.353
X0=[0,1,0]; % Lorenz's initial value
sigm=10; % Prandtl number, from Lorenz
b=8/3;   % unnamed aspect ratio, from Lorenz
h=0.01;   % Use same time step for same accuracy
steps=5000; % Complete trajectory to t=50.
t=h*(0:steps); % time coordinate

%%%% (a) Trajectories for r=166.3
rayl=166.3; 
Xdot = @(V) [sigm*(V(2)-V(1)), (rayl-V(3))*V(1)-V(2), V(1)*V(2)-b*V(3)];
X = rk4b(X0, Xdot, steps,h);  % Complete trajectory to t=50.
%%% Individual components of a trajectory on new graphs
figure; plot(t, X(:,1)); title("X Coordinate, r=166.3"); xlabel("Time"); 
figure; plot(t, X(:,2)); title("Y Coordinate, r=166.3"); xlabel("Time"); 
figure; plot(X(:,1),X(:,2)); title("(X,Z) Projection, r=166.3"); xlabel("X"); ylabel("Z"); 

%%%% (b) Trajectories for r=212
rayl=212; 
Xdot = @(V) [sigm*(V(2)-V(1)), (rayl-V(3))*V(1)-V(2), V(1)*V(2)-b*V(3)];
X = rk4b(X0, Xdot, steps,h);  % Complete trajectory to t=50.
%%% Individual components of a trajectory on new graphs
figure; plot(t, X(:,1)); title("X Coordinate, r=212"); xlabel("Time"); 
figure; plot(t, X(:,2)); title("Y Coordinate, r=212"); xlabel("Time"); 
figure; plot(X(:,1),X(:,2)); title("(X,Z) Projection, r=212"); xlabel("X"); ylabel("Z"); 

%%%% (c) Trajectories for r=146, 150, 165
rayl=146; 
Xdot = @(V) [sigm*(V(2)-V(1)), (rayl-V(3))*V(1)-V(2), V(1)*V(2)-b*V(3)];
X = rk4b(X0, Xdot, steps,h);  % Complete trajectory to t=50.
%%% Individual components of a trajectory on new graphs
figure; plot(t, X(:,1)); title("X Coordinate, r=146"); xlabel("Time"); 
figure; plot(t, X(:,2)); title("Y Coordinate, r=146"); xlabel("Time"); 
figure; plot(X(:,1),X(:,2)); title("(X,Z) Projection, r=146"); xlabel("X"); ylabel("Z"); 

rayl=151; 
Xdot = @(V) [sigm*(V(2)-V(1)), (rayl-V(3))*V(1)-V(2), V(1)*V(2)-b*V(3)];
X = rk4b(X0, Xdot, steps,h);  % Complete trajectory to t=50.
%%% Individual components of a trajectory on new graphs
figure; plot(t, X(:,1)); title("X Coordinate, r=151"); xlabel("Time"); 
figure; plot(t, X(:,2)); title("Y Coordinate, r=151"); xlabel("Time"); 
figure; plot(X(:,1),X(:,2)); title("(X,Z) Projection, r=151"); xlabel("X"); ylabel("Z"); 

rayl=165; 
Xdot = @(V) [sigm*(V(2)-V(1)), (rayl-V(3))*V(1)-V(2), V(1)*V(2)-b*V(3)];
X = rk4b(X0, Xdot, steps,h);  % Complete trajectory to t=50.
%%% Individual components of a trajectory on new graphs
figure; plot(t, X(:,1)); title("X Coordinate, r=165"); xlabel("Time"); 
figure; plot(t, X(:,2)); title("Y Coordinate, r=165"); xlabel("Time"); 
figure; plot(X(:,1),X(:,2)); title("(X,Z) Projection, r=165"); xlabel("X"); ylabel("Z"); 


%%% Ex.7
%  Function to compute Cobweb graph coordinates
function [X,Y] = cobweb(F,X0,steps)
% Input
%  F  = function handle, computes X(n+1) from X(n)
%  X0 = initial value X(1)
%  steps = number of visits to X=Y line (e.g., 10)
% Output
%  X,Y = coordinates to plot a cobweb starting at X0
  X(1) = X0;  % Initial value for iteration
  Y(1) = 0;   % Initial cobweb point is on X-axis
  for j=2:2:2*steps  % roundtrip to X=Y takes 2 steps
    X(j) = X(j-1);   % up/down to the graph of F
    Y(j) = F(X(j));  % 
    X(j+1) = Y(j);   % left/right to X=Y line
    Y(j+1) = Y(j);   % 
  end
end
%% Iteration function f(x) = 3x-x^3;
F = @(x) 3*x-x.*x.*x; 
%% (a) Cobweb graph from x0=1.9;
[X,Y]=cobweb(F, 1.9, 30);
figure;                  % New plot
t=min(X):0.01:max(X); plot(t,t);   % line (x,y)=(t,t)
hold on;  plot(t,F(t));  % iteration function
plot(X, Y);
title("Cobweb for Iteration of 3x-x^3 from x0=1.9");
xlabel("X(n)"); ylabel("X(n+1)"); 
%% (b) Cobweb graph from x0=2.1;
[X,Y]=cobweb(F, 2.1, 2);  % 2 steps
figure;                  % New plot
t=min(X):0.01:max(X); plot(t,t);   % line (x,y)=(t,t)
hold on;  % plot(t,F(t));  % iteration function omitted
plot(X, Y);
title("Cobweb for Iteration of 3x-x^3 from x0=2.1");
xlabel("X(n)"); ylabel("X(n+1)"); 

[X,Y]=cobweb(F, 2.1, 3);  % 3 steps
figure;                  % New plot
t=min(X):0.01:max(X); plot(t,t);   % line (x,y)=(t,t)
hold on;  % plot(t,F(t));  % iteration function omitted
plot(X, Y);
title("Cobweb for Iteration of 3x-x^3 from x0=2.1");
xlabel("X(n)"); ylabel("X(n+1)"); 

%%% Ex.9
%  Function to compute Cobweb graph coordinates
function [X,Y] = cobweb(F,X0,steps)
% Input
%  F  = function handle, computes X(n+1) from X(n)
%  X0 = initial value X(1)
%  steps = number of visits to X=Y line (e.g., 10)
% Output
%  X,Y = coordinates to plot a cobweb starting at X0
  X(1) = X0;  % Initial value for iteration
  Y(1) = 0;   % Initial cobweb point is on X-axis
  for j=2:2:2*steps  % roundtrip to X=Y takes 2 steps
    X(j) = X(j-1);   % up/down to the graph of F
    Y(j) = F(X(j));  % 
    X(j+1) = Y(j);   % left/right to X=Y line
    Y(j+1) = Y(j);   % 
  end
end
%% Iteration function f(x) = x(1-x), Logistic with r=1;
F = @(x) x.*(1-x);   % use componentwise ".*", for graphing
%% (a) Cobweb graph from x0=0.1;
[X,Y]=cobweb(F, 0.1, 30);
figure;                  % New plot
t=min(X):0.01:max(X); plot(t,t);   % line (x,y)=(t,t)
hold on;  plot(t,F(t));  % iteration function
plot(X, Y);
title("Cobweb for Iteration of x(1-x) from x0=0.1");
xlabel("X(n)"); ylabel("X(n+1)");
%% (b) Cobweb graph from x0=0.2;
[X,Y]=cobweb(F, 0.2, 30);
figure;                  % New plot
t=min(X):0.01:max(X); plot(t,t);   % line (x,y)=(t,t)
hold on;  plot(t,F(t));  % iteration function
plot(X, Y);
title("Cobweb for Iteration of x(1-x) from x0=0.2");
xlabel("X(n)"); ylabel("X(n+1)");
%% (c) Cobweb graph from x0=0.5;
[X,Y]=cobweb(F, 0.5, 30);
figure;                  % New plot
t=min(X):0.01:max(X); plot(t,t);   % line (x,y)=(t,t)
hold on;  plot(t,F(t));  % iteration function
plot(X, Y);
title("Cobweb for Iteration of x(1-x) from x0=0.5");
xlabel("X(n)"); ylabel("X(n+1)");
%% (d) Cobweb graph from x0=0.9;
[X,Y]=cobweb(F, 0.9, 30);
figure;                  % New plot
t=min(X):0.01:max(X); plot(t,t);   % line (x,y)=(t,t)
hold on;  plot(t,F(t));  % iteration function
plot(X, Y);
title("Cobweb for Iteration of x(1-x) from x0=0.9");
xlabel("X(n)"); ylabel("X(n+1)");
%% (e) Cobweb graph from x0=0.03;
[X,Y]=cobweb(F, 0.03, 30);
figure;                  % New plot
t=min(X):0.01:max(X); plot(t,t);   % line (x,y)=(t,t)
hold on;  plot(t,F(t));  % iteration function
plot(X, Y);
title("Cobweb for Iteration of x(1-x) from x0=0.03");
xlabel("X(n)"); ylabel("X(n+1)");
%% (f) Cobweb graph from x0=0.01;
[X,Y]=cobweb(F, 0.01, 30);
figure;                  % New plot
t=min(X):0.01:max(X); plot(t,t);   % line (x,y)=(t,t)
hold on;  plot(t,F(t));  % iteration function
plot(X, Y);
title("Cobweb for Iteration of x(1-x) from x0=0.01");
xlabel("X(n)"); ylabel("X(n+1)");


%%% Ex.10
%  Compute and plot the orbit diagram for logistic maps. 
%  (3) Best and final plot:
rmin = 3.4;  rmax = 3.7;  % range of parameter values
dr = 0.0001;               % parameter increment
rs = rmin:dr:rmax;
x = 0.33;                % initial value for iteration
ntrans = 1000;           % iterations to pass transients
norbit = 300;            % iterations within an orbit
n=0;        % count the number of points plotted
for(j=1:length(rs))   % increment over parameter range
 % Iterate with a given r, always starting with last x
 r=rs(j);
 for(i=1:ntrans)  % transient portion
   x=r*x*(1-x);
 end
 for(i=1:norbit)   % orbit portion
   n = n+1;
   x = r*x*(1-x);
   xo(n)=x;     % orbit point
   ro(n)=r;     % r value for this orbit point
   end
end
figure; hold on;   % New graph, accumulate points
plot(rmin,0,"."); plot(rmax,1,".");  % plot limits
plot(ro,xo,".");         % plot the orbit points
title("Orbit Diagram for Logistic Map Iteration");
xlabel("r, in r*x*(1-x)"); ylabel("Orbit points"); 
%%%%%  Earlier versions: %%%%%
% (1) Fewer r samples and fewer orbit points
clear("ro","xo","rs");   % delete the previous arrays
dr = 0.002;               % parameter increment  NEW!
rs = rmin:dr:rmax;
x = 0.33;                % initial value for iteration
ntrans = 1000;           % iterations to pass transients
norbit = 32;            % iterations within an orbit  NEW!
n=0;        % count the number of points plotted
for(j=1:length(rs))   % increment over parameter range
 % Iterate with a given r, always starting with last x
 r=rs(j);
 for(i=1:ntrans)  % transient portion
   x=r*x*(1-x);
 end
 for(i=1:norbit)   % orbit portion
   n = n+1;
   x = r*x*(1-x);
   xo(n)=x;     % orbit point
   ro(n)=r;     % r value for this orbit point
   end
end
figure; hold on;   % New graph, accumulate points
plot(rmin,0,"."); plot(rmax,1,".");  % plot limits
plot(ro,xo,".");         % plot the orbit points
title("Orbit Diagram for Logistic Map Iteration");
xlabel("r, in r*x*(1-x)"); ylabel("Orbit points"); 
% (2) Fewer r samples but more orbit points
clear("ro","xo","rs");   % delete the previous arrays
dr = 0.002;               % parameter increment  NEW!
rs = rmin:dr:rmax;
x = 0.33;                % initial value for iteration
ntrans = 1000;           % iterations to pass transients
norbit = 300;            % iterations within an orbit NEW!
n=0;        % count the number of points plotted
for(j=1:length(rs))   % increment over parameter range
 % Iterate with a given r, always starting with last x
 r=rs(j);
 for(i=1:ntrans)  % transient portion
   x=r*x*(1-x);
 end
 for(i=1:norbit)   % orbit portion
   n = n+1;
   x = r*x*(1-x);
   xo(n)=x;     % orbit point
   ro(n)=r;     % r value for this orbit point
   end
end
figure; hold on;   % New graph, accumulate points
plot(rmin,0,"."); plot(rmax,1,".");  % plot limits
plot(ro,xo,".");         % plot the orbit points
title("Orbit Diagram for Logistic Map Iteration");
xlabel("r, in r*x*(1-x)"); ylabel("Orbit points");