Crank-Nicolson

function [error,h,k] = heateq_CN(a,b,ti,tf,m,kappa)
% a = 0;
% b = 1;
% m = 100;
% ti = 0;%initial time
% tf = 1.1;%final_time
% kappa = 0.02;
    % Discretization set up
    h = abs(b-a)/(m+1);            % h = delta x               
    k = 4*h;                       % time step k
    x = linspace(a,b,m+2)';        % x(1)= a and x(m+2)= b    
    step = round((tf-ti)/k);       % approx. number of time steps

    % Analytical solution 
    u_true = @(x,t) 1 + sin(pi*x).*exp(-0.02*pi^2*t) + sin(3*pi*x).*exp(-0.02*9*pi^2*t);
    u0 = u_true(x,0);              % Initial condition
 
    % Matrix Form
    r = (1/2) * kappa* k/(h^2);    % Define the parameter
    e = ones(m,1);
    S = spdiags([e -2*e e], [-1 0 1], m, m);
    A = eye(m) - r * S;            % Tri- diagonal matrix
    A1 = eye(m) + r * S;           % Matrix for rhs
    
    t_n = 0; u = u0; 
    U_soln = []; t = 0;
    U_soln(:,1) = u;

    % Time loops
    for i = 1 : step        
         t_np = t_n + k;   % t_{n+1}

         % Take the interior data
         u_in = u;         
         u_in(1) = [];
         u_in(end) = [];
         
         % Calculate the rhs vector
         rhs = A1*u_in;    
         rhs(1) = rhs(1) + 2*r;
         rhs(m) = rhs(m) + 2*r;
         
         % Solve the problem at time t_{n+1}
         u_in = A\rhs; 
         u = [1; u_in; 1];   % Add boundary profile
         U_soln(:,i+1) = u;  % Record Solution
         t_n = t_np;   
         t(end+1) = t_n;     % Record time 
    end
         error = max(abs(u-u_true(x,t_np)));
         clf
         [X,Y] = meshgrid(x,t);
         mesh(X,Y,U_soln');   % Plot the temperature distribution 
         xlabel('x');ylabel('t');
         title('Numerical solution to the heat flow problem when h = 0.1, k = 0.4');  
         
end