function homework1_solutions

% MAX TOTAL GRADE 30 POINTS 

close all   % This closes any open fig windows

% Solution to Problem #1   [2 POINTS]
for i=1:401
   x(i) = -6+12*(i-1)/400;
end
% Or, use x = linspace(-6,6,401)
% Or      x = [-6:6:12/400]

% Solution to problem #2   [2 POINTS]
figure1 = figure;
for i=1:length(x)
  y(i) = sin(8*x(i))*exp(-x(i)^2);
end
% or use y = sin(8*x).*exp(-x.^2)
% The axis labels etc are optional for all problems...
axes('Parent',figure1,'FontSize',12);
box('on');
hold('all');
plot(x,y);
title({'Solution to problem #2'},'FontSize',14);
xlabel({'{\it x}'},'FontSize',12,'FontName','Times New Roman');
ylabel({'y'},'FontSize',12,'FontName','Times New Roman',...
    'FontAngle','italic');


% Solution to problem #3  [3 POINTS]
for i = 1:101
  v(i) = (i-1)*2*pi/100;
  sinv = sin(v(i));
  if (sinv<-0.707)
      w(i) = -0.7071;
  elseif (sinv>0.707)
      w(i) = 0.707;
  else
      w(i) = sinv;
  end
end
figure2 = figure;
axes('Parent',figure2,'FontSize',12);
box('on');
hold('all');
plot(v,w);
title({'Solution to problem #3'},'FontSize',14);
xlabel({'{\it v}'},'FontSize',14,'FontName','Times New Roman');
ylabel({'w'},'FontSize',14,'FontName','Times New Roman',...
    'FontAngle','italic');

% Solution to problem #4 (see below for the function) [5 POINTS]
rms = compute_rms(w)
%The numbers you get for the rms depend on the number of
% points

%Solution to Problem #5 [5 POINTS]
beta = 0.8;   % Maximal expression rate, mols/(sec*m^3)
alpha =0.6;   % Degredation rate, 1/sec.
K = 0.25;   % Activation coefficient, mols/m^3
n = 4;
C0 = 0.1;   % Initial concentration (mols/m^3)
[tvals1,cvals1] = ode45(@protein,[0,15.],C0);
C0 = 0.2;   % Initial concentration (mols/m^3)
[tvals2,cvals2] = ode45(@protein,[0,15.],C0);
figure3 = figure;
axes('Parent',figure3,'FontSize',12);
box('on');
hold('all');
plot(tvals1,cvals1);
plot(tvals2,cvals2);
title({'Solution to problem #5'},'FontSize',14);
xlabel({'{Time (sec)}'},'FontSize',14,'FontName','Times New Roman');
ylabel({'% Concentration (mols/m^3)'},'FontSize',14,'FontName','Times New Roman');

% 
% When the concentration reaches steady-state, dC/dt=0.
% The number of possible solutions depends on the number of
% roots of the equation
%  beta*(C/K)^n/(1+(c/K)^n)   - alpha*C=0
% If n=1, then the equation reduces to 
% beta*(C/K)/(1+(c/K))   - alpha*C=0
% => C*((beta/K-alpha)   - alpha*C/K)=0
%  This has solutions C=0 (which is unstable for the values of 
%  beta, K and alpha used in the problem - a small change in the
%  concentration will make it increase to the stable point.  You can
%  show this with a computation -e.g. start with C=10^-9 )
%  and C= (beta-alpha*K)/alpha  (which is stable)
%  [5 POINTS]


% Solution to problem #6 [8 POINTS]

% Parameters for 6.1
x0 = -0.2706;  
y0 = -3.3019; 
z0 = 70.3189; 
rho=99.96;
w0 = [x0,y0,z0];
[timevals1,solvals1] = ode45(@lorenz,[0,20],w0);
figure4 = figure;
axes61 = axes('Parent',figure4,'FontSize',12);
box('on');
hold('all');
plot61 = plot(timevals1,solvals1);
set(plot61(1),'DisplayName','x');
set(plot61(2),'DisplayName','y');
set(plot61(3),'DisplayName','z');
title({'Solution to problem #6.1'},'FontSize',14);
xlabel({'{Time (sec)}'},'FontSize',14,'FontName','Times New Roman');
ylabel({'Solution'},'FontSize',14,'FontName','Times New Roman');
legend61 = legend(axes61,'show');
set(legend61,'Position',[0.75 0.7492 0.1268 0.1587]);

% Problem 6.2
% Parameters for 6.1
x0 = -0.2706;  
y0 = -3.3019; 
z0 = 70.3189; 
rho=28;
w0 = [x0,y0,z0];
[timevals2,solvals2] = ode45(@lorenz,[0,20],w0);

figure5 = figure;
axes62 = axes('Parent',figure5,'FontSize',12);
box('on');
hold('all');
plot62 = plot(timevals2,solvals2);
set(plot62(1),'DisplayName','x');
set(plot62(2),'DisplayName','y');
set(plot62(3),'DisplayName','z');
title({'Solution to problem #6.2'},'FontSize',14);
xlabel({'{Time (sec)}'},'FontSize',14,'FontName','Times New Roman');
ylabel({'Solution'},'FontSize',14,'FontName','Times New Roman');
legend62 = legend(axes62,'show');
set(legend62,'Position',[0.75 0.7492 0.1268 0.1587]);

figure6 = figure;
plot3(solvals1(:,1),solvals1(:,2),solvals1(:,3))
xlabel({'x'});
ylabel({'y'});
zlabel({'z'});

figure7 = figure;
plot3(solvals2(:,1),solvals2(:,2),solvals2(:,3))

xlabel({'x'});
ylabel({'y'});
zlabel({'z'});


    function dCdt = protein(t,C)
        dCdt = beta*(C/K)^n/(1+(C/K)^n) - alpha*C;
    end

    function dwdt = lorenz(t,w)
        % The array w stores [x,y,z] at time t
        x = w(1); y = w(2); z=w(3);
        dxdt = 10*(y-x);
        dydt = x*(rho-z)-y;
        dzdt = x*y-8*z/3;
        dwdt = [dxdt;dydt;dzdt];
    end

end
function rms = compute_rms(vector)
% Function to compute the rms value of a vector
  rms = 0;
  for i =1:length(vector)
      rms = rms + vector(i)*vector(i);
  end
  rms = sqrt(rms/length(vector));
  % Or, faster but more obscure rms = norm(vector)/sqrt(length(vector));
end