MATLAB: Problem with Index (Dufort Frankel Scheme)

correct formatdufort frankelindex

I am currently trying to code this to fit a Dufort Frankel Scheme. The scheme is attached below as an image. This is the code I have currenlty and I am getting several errors, and I am not sure this is correct. Can someone help with this problem? The index is exceeding the array bounds in line 55 and is the correct way to write the Dufort Scheme out? 
clear all;                  % clear all variables in memory
 
xl=0; xr=1;                 % x domain [xl,xr]
J = 10;                     % J: number of division for x
dx = (xr-xl) / J;           % dx: mesh size
tf = 0.1;                   % final simulation time
Nt = 50;                    % Nt: number of time steps
dt = tf/Nt; 
k = 0;
t = 0;
 
mu = dt/(dx)^2; 
 
if mu > 0.5         % make sure dt satisy stability condition
    error('mu should < 0.5!')
end
% Evaluate the initial conditions
x = xl : dx : xr;              % generate the grid point
% f(1:J+1) since array index starts from 1
f = sin(pi*x) + sin(2*pi*x);  
 
% store the solution at all grid points for all time steps
u = zeros(J+1,Nt);   
 
% Find the approximate solution at each time step
for n = 1:Nt
    t = n*dt;         % current time
   
    % boundary condition at left side
    gl = sin(pi*xl)*exp(-pi*pi*t)+sin(2*pi*xl)*exp(-4*pi*pi*t);   
    % boundary condition at right side
    gr = sin(pi*xr)*exp(-pi*pi*t)+sin(2*pi*xr)*exp(-4*pi*pi*t);  
    if n==1    % first time step
       for j=2:J    % interior nodes     
       u(j,n) = f(j) + mu*2*(f(j+1)-(f(j)+f(j))+f(j-1));
       end     
       u(1,n) = gl;   % the left-end point

       u(J+1,n) = gr; % the right-end point 

    else 
       for j=2:J    % interior nodes
         u(j,n)=u(j,n-1)+mu*2*(u(j+1,n)-(u(j,n+1)+u(j,n-1))+u(j-1,n));
       end
       u(1,n) = gl;   % the left-end point
       u(J+1,n) = gr; % the right-end point 
    end
    
    % calculate the analytic solution 
   
    for j=1:J+1
        xj = xl + (j-1)*dx;
        u_ex(j,n)=sin(pi*xj)*exp(-pi*pi*t) ...
+sin(2*pi*xj)*exp(-4*pi*pi*t);
    end
end
 
% Plot the results
tt = dt : dt : Nt*dt;
figure(1)
colormap(gray);     % draw gray figure
surf(x,tt, u');     % 3-D surface plot

xlabel('x')
ylabel('t')
zlabel('u')
title('Numerical solution of 1-D parabolic equation')
 
figure(2)
surf(x,tt, u_ex');     % 3-D surface plot
xlabel('x')
ylabel('t')
zlabel('u')
title('Analytic solution of 1-D parabolic equation')
 
maxerr=max(max(abs(u-u_ex)));

Best Answer

I strongly advice you to check and rethink on your code. 
xl=0; xr=1;                 % x domain [xl,xr]
J = 10;                     % J: number of division for x
dx = (xr-xl) / J;           % dx: mesh size
tf = 0.1;                   % final simulation time
Nt = 50;                    % Nt: number of time steps
dt = tf/Nt;
k = 0;
t = 0;
mu = dt/(dx)^2;
if mu > 0.5         % make sure dt satisy stability condition
    error('mu should < 0.5!')
end
% Evaluate the initial conditions
x = xl : dx : xr;              % generate the grid point
% f(1:J+1) since array index starts from 1
f = sin(pi*x) + sin(2*pi*x);
% store the solution at all grid points for all time steps
u = zeros(J+1,Nt);
u_ex = zeros(J+1,Nt) ;
% Find the approximate solution at each time step
for n = 1:Nt-1
    t = n*dt;         % current time
    
    % boundary condition at left side
    gl = sin(pi*xl)*exp(-pi*pi*t)+sin(2*pi*xl)*exp(-4*pi*pi*t);
    % boundary condition at right side
    gr = sin(pi*xr)*exp(-pi*pi*t)+sin(2*pi*xr)*exp(-4*pi*pi*t);
    if n==1    % first time step
        for j=2:J    % interior nodes

            u(j,n) = f(j) + mu*2*(f(j+1)-(f(j)+f(j))+f(j-1));
        end
        u(1,n) = gl;   % the left-end point

        u(J+1,n) = gr; % the right-end point

    else
        for j=2:J    % interior nodes
            u(j,n)=u(j,n-1)+mu*2*(u(j+1,n)-(u(j,n+1)+u(j,n-1))+u(j-1,n));
        end
        u(1,n) = gl;   % the left-end point
        u(J+1,n) = gr; % the right-end point
    end
    
    % calculate the analytic solution
    
    for j=1:J+1
        xj = xl + (j-1)*dx;
        u_ex(j,n)=sin(pi*xj)*exp(-pi*pi*t) ...
            +sin(2*pi*xj)*exp(-4*pi*pi*t);
    end
end
% Plot the results
tt = dt : dt : Nt*dt;
figure(1)
colormap(gray);     % draw gray figure
surf(x,tt, u');     % 3-D surface plot

xlabel('x')
ylabel('t')
zlabel('u')
title('Numerical solution of 1-D parabolic equation')
figure(2)
surf(x,tt, u_ex');     % 3-D surface plot
xlabel('x')
ylabel('t')
zlabel('u')
title('Analytic solution of 1-D parabolic equation')
maxerr=max(max(abs(u-u_ex)));
Related Question