heat equation with Neumann B.C in matlab

Question

$$ \frac{\partial u}{\partial t}=\alpha\frac{\partial^{2}u}{\partial x^{2}} \qquad u(x,0)=f(x)\qquad u_{x}(0,t)=0\qquad u_{x}(1,t)=2 $$ i'm trying to code the above heat equation with neumann b.c. using explicit forward finite differences in matlab. I used central finite differences for boundary conditions. $$ u_{x}(0,t)=\frac{u_{i+1}^{j}-u_{i-1}^{j}}{2h} $$ for i=1 ı used $$ u_{2}^{j}-u_{x}(0,t)2h= u_{0}^{j} $$ and for i=m $$ u_{m-1}^{j}-u_{x}(0,t)2h= u_{m+1}^{j} $$ Here my code and it doesn't give correct results. Actually i am not sure that i coded correctly the boundary conditions. I call the function as heatNeumann(0,0.1,0,1,0.1,0.0001,1) It would be good if someone can help.

function heatNeumann(t_i,t_f,a,b,dx,dt,alpha)

% dt: step size in time dimension
% dx: step size in x axes
% a: left point of domain x
% b: right point of domain x
% alpha: take 1
% t_i: t initial, t_f: t final



n=(t_f-t_i)/dt;
m=(b-a)/dx;
lambda=alpha*dt/dx^2;
T=zeros(m+1,n+1);



x=a:dx:b;
t=t_i:dt:t_f;

u0 = x.^2+1+cos(pi.*x);
T(:,1)=u0; %initial value
u_left = T(2,1); %boundary cond.
u_right = T(m,1)+4*dx; %boundary cond.

for j=1:n
    T(1,j+1)=T(1,j)+lambda*(T(2,j)-2*T(1,j)+u_left); %boundary cond.
    for i=2:m
        T(i,j+1)=T(i,j)+lambda*(T(i+1,j)-2*T(i,j)+T(i-1,j));
    end
    T(m+1,j+1)=T(m+1,j)+lambda*(u_right-2*T(m+1,j)+T(m,j)); %boundary cond.
    u_left = T(2,j+1);
    u_right = T(m,j+1);
end


%exact soln

[xx,tt]=meshgrid(x,t);
exact=2.*tt+xx.^2+1+exp(-pi^2.*tt).*cos(pi.*xx);


[T(:,end) exact(end,:)']

Answer 1

You have the right idea, your boundary condition is, $$u_x(t_n,x_0)=\frac{v_1^m-v_{-1}^m}{2h}$$ Now apply your scheme to get $v_0^{m+1}$. That is, $$v_0^{m+1} = v_0^m + b\mu[v_1^m - 2v_0^m + v_{-1}^m]=v_0^m + b\mu\left[v_1^m - 2v_0^m + \left(v_1^m-2hu_x\left(t_n,x_0\right)\right)\right]$$

And do the same for the right boundary condition. Then your BCs should become,

T(1,j+1)=T(1,j)+lambda*(T(2,j)-2*T(1,j)+T(2,j));
T(m+1,j+1)=T(m+1,j)+lambda*(T(m,j)-2*T(m+1,j)+(T(m,j)+2*dx*2));