clear all

  N_MC = 1e+4;
  time = 12;
  dt = .01;
  tvals = dt:dt:time;
  Nt = length(tvals);
  alpha = 1;
  beta = 10;
  sigma_a = sqrt(1/48);
  sigma_b = sqrt(3);
  a1 = 3/4;
  a2 = 5/4;
  b1 = 7;
  b2 = 13;

  mean_true = beta*(sinh(sqrt(3)*sigma_a*tvals)./(sqrt(3)*sigma_a*tvals)).*exp(-alpha*tvals);
  factor1 = (beta^2 + sigma_b^2)*sinh(2*sqrt(3)*sigma_a*tvals)./(2*sqrt(3)*sigma_a*tvals);
  factor2 = - beta^2* (sinh(sqrt(3)*sigma_a*tvals).^2)./(3*sigma_a^2*tvals.^2);
  var_true = (factor1 + factor2).*exp(-2*alpha*tvals);
  mean_true = (40./tvals).*sinh(tvals/4).*exp(-tvals); 
  var_true = ((206./tvals).*sinh(tvals/2) - (1600./tvals.^2).*(sinh(tvals/4).^2) ).*exp(-2*tvals);
  sd_true = sqrt(var_true);
  mean2p_true = mean_true + 2*sd_true;
  mean2m_true = mean_true - 2*sd_true;

  figure(1)
  plot(tvals,mean_true,'-b',tvals,mean2p_true,'-b',tvals,mean2m_true,'-b','linewidth',2)
  %plot(tvals,mean_true,'-b',tvals,mean2p_true,'-b',tvals,mean2m_true,'-b',tvals,0*tvals,'linewidth',2)
  set(gca,'Fontsize',[20]);
  xlabel('Time (s)')
  ylabel('Displacement (cm)')
%  legend('True Mean','True 2\sigma Interval','Location','NorthEast')


%
% Compute N_MC samples for alpha and beta
%

  alpha_MC = a1 + (a2-a1)*rand(N_MC,1);
  beta_MC = b1 + (b2-b1)*rand(N_MC,1);  

  for k=1:Nt
    MC = beta_MC.*exp(-alpha_MC*tvals(k));
    MC_mean(k) = (1/N_MC)*sum(MC);
    tmp = 0;
    for j=1:N_MC
      tmp = tmp + (MC(j) - MC_mean(k))^2;
    end
    MC_sigma(k) = sqrt((1/(N_MC-1))*tmp);
  end  

  figure(2)
  plot(tvals,mean_true,'-b',tvals,MC_mean,'--k');

  figure(3)
  plot(tvals,sd_true,'-b',tvals,MC_sigma,'--k');