As I can put the Neumann boundary conditions? the Crank Nicholson scheme, the Heat Equation

Question

How i can put the Neumann BC in my code? I tested but I get error, because the arrays are not the same size

$$U_t=U_{xx},\quad 0<x<1$$

$$u_x(0,t)=0$$ $$u_x(1,t)=0$$

$$u(x,0)=f(x)$$

I have my equation discretized by the method and my program

the BC Neumann it are when $i=1$ :

$$u(2,j)=u(0,j)$$

$$u(n+1,j)=u(n-1,j)$$

here is the problem with the boths BC

an apology if not in code here is somewhat long, thanks for your help!

Answer 1

By my comment, here is what I meant:

I solved the PDE analytically and obtained: $$ u(x, t) = \frac{1}{\pi}\sum_{n = 1}^{\infty}\frac{\sin[(2n -1)\pi x]}{2n -1}\exp(-\pi^2(2n-1)^2t) $$ where I assumed $f(x) = 1$ just so I could plot something.

% PDE

clear all
close all

Nmax = 500;                            % truncate series
x = linspace(0, 1, 100000);            % x vector
t = [0, 0.0001, 0.001, 0.01, 1];       % time steps to plot
U = zeros(length(x), 1)';              % pre-allocate U

cc = hsv(length(t));                   % different colors per line

for k = 1:length(t)
    for n = 1:Nmax
        U = U + 1/(pi*(2*n - 1))*sin((2*n - 1)*pi*x)...
            *exp(-pi^2*(2*n-1)^2*t(k));
    end
    hold on
    plot(x, U, 'color', cc(k,:))
end