clear all
  
  N = input('N = ');
  beta = 10;
  mu = 1.0;
  sigma = 0.25;
  sigma_beta = 2;
  sigma_beta = 0;
  time = 12;
  dt = .01;
  Ntime = time/dt;

  b = zeros(N+1,1);
  b(1,1) = beta;
  b(2,1) = sigma_beta;

  j_ind1 = 0;
  j_ind2 = 1;
  for i = 1:N+1
    i_ind = i-1;
    h_i = factorial(i_ind); 
    for k = 1:N+1
      k_ind = k-1;
      s2_ind1 = i_ind + j_ind1 + k_ind;
      s2_ind2 = i_ind + j_ind2 + k_ind;
      if (mod(s2_ind1,2) == 0) & (s2_ind1/2 >= i_ind) & ( s2_ind1/2 >= k_ind) & (s2_ind1/2 >= j_ind1)
        s1 = s2_ind1/2;
        num1 = factorial(i_ind)*factorial(j_ind1)*factorial(k_ind);
        den1 = factorial(s1-i_ind)*factorial(s1-j_ind1)*factorial(s1-k_ind);
        e1 = num1/den1;
      else
        e1 = 0;
      end
      if (mod(s2_ind2,2) == 0) & (s2_ind2/2 >= i_ind) & ( s2_ind2/2 >= k_ind) & (s2_ind2/2 >= j_ind2)
        s2 = s2_ind2/2;
        num2 = factorial(i_ind)*factorial(j_ind2)*factorial(k_ind);
        den2 = factorial(s2-i_ind)*factorial(s2-j_ind2)*factorial(s2-k_ind);
        e2 = num2/den2;
      else
        e2 = 0;
      end
      A(i,k) = -(1/h_i)*(mu*e1 + sigma*e2);
    end
  end

  AA = inv(eye(N+1,N+1) - dt*A);
  U = b;
  tvals(1) = 0;
  
  var(1) = sigma_beta^2;
  for n = 1:Ntime
    tvals(n+1) = n*dt; 
    U(:,n+1) = AA*U(:,n);
%    U(:,n+1) = (eye(N+1,N+1) + dt*A)*U(:,n);
    var_sum = 0;
    for j=2:N+1
      var_sum = var_sum + (U(j,n+1)^2)*factorial(j-1);
    end 
    var(n+1) = var_sum;
  end

  mean = U(1,:);   
  sd = sqrt(var);
  mean2p = mean + 2*sd;
  mean2m = mean - 2*sd;
  det = exp(-mu*tvals)*beta; 
  mean_true = beta*exp(-mu*tvals).*exp(.5*(sigma^2)*(tvals.^2));
  var_true = exp(-2*mu*tvals).*(exp(2*(sigma^2)*(tvals.^2))*(beta^2 + sigma_beta^2) - exp((sigma^2)*(tvals.^2))*beta^2);
  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,mean,'--b',tvals,mean2p_true,'-b',tvals,mean2p,'--b', ...
       tvals,mean2m_true,'-b',tvals,mean2m,'--b','linewidth',2)
  set(gca,'Fontsize',[20]);
  xlabel('Time (s)')
  ylabel('Displacement (cm)')
  legend(' Analytic',' Approximate','Location','NorthEast')