function tdiffusion1d_rk4_vs_samarskii (k, tTime, max_time_steps, nElem)

  clc; format short g;
  
  %%% Pure Diffusion Problem (Transient)
  %%% Study of Runge--Kutta explicit method
  %%% From the desk of Prashanth Nadukandi
  
  %%% Problem data (1D transient diffusion)
  L = 1.0; % 1d domain: [-L, L]
  % k = 1e-3;
  blink = 1e-2; % pause between subsequent solution plot
  
  %%% Discretization data
  % nElem           = 334; % trouble starts at a mesh size > 333
  % tTime           = 5;
  % max_time_steps  = 200;
  
  nNode = nElem+1;
  nSize = nNode-2; % as the two boundary node values will be prescribed.
  dL = 2*L/nElem;
  dT = tTime/max_time_steps;
  x_h = linspace(-L,L,nNode);
  x_a = linspace(-L,L,2001);
  
  % System matrices
  I = speye(nSize, nSize); % sparse identity matrix
  A = gallery('tridiag',ones(nSize-1,1),-2*ones(nSize,1),ones(nSize-1,1)); % LHS
  A = (k/(dL^2))*A;
  b = zeros(nSize,1); % RHS; boundary conditions will be included later
  
  % Apply the Samarskii's additive operator splitting (AOS) on A
  % B := tril(A,-1) + diag(diag(A))/2
  % C := triu(A, 1) + diag(diag(A))/2
  B = gallery('tridiag',ones(nSize-1,1),-ones(nSize,1),0*ones(nSize-1,1));
  C = gallery('tridiag',0*ones(nSize-1,1),-ones(nSize,1),ones(nSize-1,1));
  B = (k/(dL^2))*B;
  C = (k/(dL^2))*C;
  
  % The Initial Solution Profile: Double Square Pulse
  Phi_0 = zeros(nNode,1);
  aa = 0.2;
  for i = 1:nNode
      % Square Pulse (to model a shock regime)
      if (x_h(i)>=-aa) && (x_h(i)<=aa)
          Phi_0(i) = 1.0;
      end
  end
  
  Phi_rk4 = Phi_0;
  Phi_sam = Phi_0;
  %%% Include the boundary conditions
  b(1) = (k/(dL^2))*Phi_0(1); b(end) = (k/(dL^2))*Phi_0(end);
  
  % Preface for plotting solution
  YMIN = min(Phi_0); YMAX = max(Phi_0);
  fig1 = cfigure(35,15);
  figure(fig1);
  %grid on;
  subplot(1,2,1), hold on;
  subplot(1,2,1), axis tight;
  subplot(1,2,1), set(gca,'XLim',[-L L]);
  subplot(1,2,1), set(gca,'YLim',[-0.1+YMIN 0.1+YMAX]);
  subplot(1,2,1), set(gca,'XTick',-L:0.2:L);
  subplot(1,2,1), set(gca,'YTick',-0.5+YMIN:0.1:0.5+YMAX);
  subplot(1,2,1), xlabel('Space');
  subplot(1,2,1), ylabel('\Phi');
  subplot(1,2,1), title('Runge-Kutta 4','FontSize',18);
  subplot(1,2,1), ...
  plot(x_h, Phi_0, '-r','Linewidth',2); % profile of the initial solution
  subplot(1,2,1), set(gcf, 'Color','white');
  
  subplot(1,2,2), hold on;
  subplot(1,2,2), axis tight;
  subplot(1,2,2), set(gca,'XLim',[-L L]);
  subplot(1,2,2), set(gca,'YLim',[-0.1+YMIN 0.1+YMAX]);
  subplot(1,2,2), set(gca,'XTick',-L:0.2:L);
  subplot(1,2,2), set(gca,'YTick',-0.5+YMIN:0.1:0.5+YMAX);
  subplot(1,2,2), xlabel('Space');
  subplot(1,2,2), ylabel('\Phi');
  subplot(1,2,2), title('Samarskii AOS','FontSize',18);
  subplot(1,2,2), ...
  plot(x_h, Phi_0, '-r','Linewidth',2); % profile of the initial solution
  subplot(1,2,2), set(gcf, 'Color','white');

  %The time integration loop
  for time_step = 1:max_time_steps
 
    % compute the RK4 solution 
    % RK4 Step 1: k_1 = 0.5*f(t_n, y_n)
    Psi_1 = 0.5*(A*Phi_rk4(2:end-1) + b); 
    % RK4 Step 2: k_2 = 0.5*f(t_n+0.5*dT, y_n + dT*k_1)
    Psi_2 = 0.5*(A*(Phi_rk4(2:end-1) + dT*Psi_1) + b); 
    % RK4 Step 3: k_3 = 0.5*f(t_n+0.5*dT, y_n + dT*k_2)
    Psi_3 = 0.5*(A*(Phi_rk4(2:end-1) + dT*Psi_2) + b); 
    % RK4 Step 4: k_4 = 0.5*f(t_n+dT, y_n + 2*dT*k_3)
    Psi_4 = 0.5*(A*(Phi_rk4(2:end-1) + 2*dT*Psi_3) + b); 
    % Solve for the solution at t_n+1
    Phi_rk4(2:end-1) = (dT/3)*(Psi_1 + 2*Psi_2 + 2*Psi_3 + Psi_4) ...
                   + Phi_rk4(2:end-1);
  
    % compute the Samarskii solution
    mu = (I - B*dT)\((I+C*dT)*Phi_sam(2:end-1) + b*dT); 
    nu = (I - C*dT)\((I+B*dT)*Phi_sam(2:end-1));
    Phi_sam(2:end-1) = (mu+nu)/2;

    % compute the analytical solution
    cTime = dT*time_step;
    Phi_a = 0.5*(erf((x_a+aa)/sqrt(4*k*cTime)) - erf((x_a-aa)/sqrt(4*k*cTime)));
    % subplot(1,2,1), 
    subplot(1,2,1), cla;
    subplot(1,2,1), plot(x_h,Phi_rk4,'-r','Linewidth',2);
    subplot(1,2,1), plot(x_a,Phi_a,'--k','Linewidth',2);
    subplot(1,2,1), pause(blink);

    subplot(1,2,2), cla;
    subplot(1,2,2), plot(x_h,Phi_sam,'-r','Linewidth',2);
    subplot(1,2,2), plot(x_a,Phi_a,'--k','Linewidth',2);
    subplot(1,2,2), pause(blink);
  
  end %End of time loop
  hold off;
  
  clear all;
end
