clear all

%
% Input the model parameters and construct time steps.
%
  dt = .1;
  k = 4.2;
  m = 1.9;
  kh = 1e-8*k;
  mh = 1e-8*m;
  kdk = k + kh;
  mdm = m + mh;
  t = 0:dt:10;
  h = 1e-8;
  tdt = t + h;
 
%
% Compute the solution and derivative.
%
  y = 2*cos(sqrt(k/m)*t);
  dydt = 2*(cos(sqrt(k/m)*tdt) - cos(sqrt(k/m)*t))/h;

%
% Compute the analytic sensitivity relations.
%

  dydmfd = 2*(cos(sqrt(k/mdm)*t) - cos(sqrt(k/m)*t))/mh;
  dydkfd = 2*(cos(sqrt(kdk/m)*t) - cos(sqrt(k/m)*t))/kh;
  dydm = sqrt(k/m^3)*t.*sin(sqrt(k/m)*t);
  dydk = -(t/sqrt(m*k)).*sin(sqrt(k/m)*t);

%
% Compute the Fisher information and compute its eigenvalues.
%

  S = [dydm' dydk'];
  V = S'*S;
  
  eig(V)

%
% Plot the analytic and finite-difference sensitivities.
%

  figure(1)
  plot(t,dydm,t,dydmfd,'--',t,0*t,'k-','linewidth',2)
  legend('Analytic','FD','Location','NorthWest')
  xlabel('Time (s)')
  ylabel('dydm')
  set(gca,'Fontsize',[20]);

  figure(2)
  plot(t,dydk,t,dydkfd,'--r',t,0*t,'k-','linewidth',2)
  legend('Analytic','FD','Location','NorthWest')
  xlabel('Time (s)')
  ylabel('dydk')
  set(gca,'Fontsize',[20]);