# Why is my BTCS diffusion PDE solver matlab code not working?

I am trying to create a function that solves one BTCS step of the following PDE:

$\frac{\partial q}{\partial t} = k\frac{\partial^{2}q}{\partial x^{2}}$

over a domain of: $0\leq x \leq \pi$

with initial condition given by the following square wave:

$q(x,0) = \left\{ \begin{array}{ll} 1, & \frac{2\pi}{5}<x<\frac{3\pi}{5}\\ 0, & otherwise \end{array} \right.$

With boundary conditions:

$q(0,t) = 0$

$\frac{\partial q}{\partial x}(\pi,t) = 0$

and $k = 0.05$

Here is my code:

function [ qnew ] = btcs_step( q,dt,dx,k )
%INPUTS:
% q: column vector of solution values from previous timestep
% dt: timestep
% dx: x-step
% k: diffusion equation coefficient

% generate A matrix
A_size = length(q);
A = gallery('tridiag', A_size, -k/dx^2, 1/dt + (2*k)/dx^2, -k/dx^2);

% generate b matrix
%b = ones(length(q),1).*(1/dt).*q';
b = q'/dt;

%impose boundary conditions
% Dirichlet conidtions
b(1) = 0;
b(end) = b(end-1);

% Modifications to coeff matrix
A(1,1) = 1;
A(1,2) = 0;
A(end,end-1) = 0;
A(end,end) = 1;

% solve matrix
qnew = A\b;
end


The second Dirichlet condition is a finite difference approximation to the derivative

I am calling this function and plotting results for various timesteps at various times with the following code:

% parameters
k = 0.05;
N = 100;

dx = pi/N;
timesteps = [0.005, 0.5*dx^2/k, 0.1];
x_vals = 0:dx:pi;

%initial condiiton
q_init = (2*pi/5)<x_vals & x_vals<(3*pi/5);

%stuff for plotting
plot_times = [0,2.5,5,10,20,40,80];
legend_text = strcat({'t = '},string(plot_times));

for dt = timesteps
% create data structure
tvals = 0:dt:80; nt = length(tvals);
q = zeros(nt, N+1);
q(1,:) = q_init;
for t = 2:nt
q(t,:) = btcs_step( q(t-1,:), dt, dx, k);
end
figure;
hold on;
for pt = plot_times
pti = floor(pt/dt)+1;
plot(x_vals,q(pti,:));
end
title(sprintf('FTCS Solution of 1-D diffusion equation\nk = 0.05\nTimestep: %g',dt));
legend(legend_text);
xticks([0,pi/4,pi/2,3*pi/4,pi]);
xticklabels({'0','\pi/4','\pi/2','3\pi/4','\pi'});
axis([0,pi,-inf,inf])
hold off
end


Here are the three plots:

There should be 7 curves on each plot...

• There are possibilities of errors at a mathematical level (using the wrong equations) and at the programming level (incorrect implementation of the equations.) You should add a discussion of the equations that you're using to your question. – Brian Borchers May 4 '18 at 15:32