function [] = Heat1d(N,m,T)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Heat1d solves heat equation using Bender-schmidt, crank-nicolson and
% Dufort-frankel methods. N: no of x axis mesh intervals, m: no of time
% intervals, T: final time
% It also compares the solutions for the equation:
% u_t-u_xx=0, u(x,0)=sin(pi x), u(0,t)=0=u(1,t)

% Exact solution: u(x,t)=sin(pi x)*exp(-pi^2*t)

% Run with typing in command window: a) Heat1d(100,100,1) b) Heat1d(50,500,0.1) c) Heat1d(10,10,0.1)
% d) Heat1d(100,10,0.1)

% Dufort frankel scheme is consistent for CFL=constant
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


dx=1/N
dt=T/m
xx=0:dx:1;
tt=0:dt:T;

uu=@(x,t) sin(pi*x).*exp(-pi^2*t);
[XX,TT]=meshgrid(xx,tt);
U=uu(XX,TT);

l=dt/dx/dx;
CFL=l
e=ones(N-1,1);

%% Bender-schmidt

% comment this section for large CFL

A=spdiags([l*e (1-2*l)*e l*e],[-1 0 1],N-1,N-1);
u=zeros(m+1,N+1);
u(1,:)=sin(pi*xx);
for n=1:m
    b=u(n,2:N)';
    u(n+1,2:N)=A*b;
end

figure(1)
mesh(xx,tt,u)

xlabel('x');
ylabel('t');
title('Solution of Heat equation : Bender-Schmidt')
Error_Bender =norm(u-U,Inf)

%% Crank-Nicolson

%e=ones(N-1,1);
A=spdiags([-l*e 2*(1+l)*e -l*e]/2,[-1 0 1],N-1,N-1);
B=spdiags([l*e 2*(1-l)*e l*e]/2,[-1 0 1],N-1,N-1);
u1=zeros(m+1,N+1);
u1(1,:)=sin(pi*xx);

for n=1:m
    b=u1(n,2:N)';
    temp=B*b;
    u1(n+1,2:N)=A\temp;
end

figure(2)
mesh(xx,tt,u1)


xlabel('x');
ylabel('t');
title('Solution of Heat equation : Crank-Nicolson')
Error_Crank =norm(u1-U,Inf)

%% Dufort-Frankel

k=2*l;r=(1-k)/(1+k);s=k/(1+k);
A=[r*eye(N-1,N-1) spdiags([s*e 0*e s*e],[-1 0 1],N-1,N-1)];
u2=zeros(m+1,N+1);
u2(1:2,:)=u1(1:2,:); % first time step is calculated by Crank-Nicolson

for n=2:m
    bb=[u2(n-1,2:N) u2(n,2:N)]';
    u2(n+1,2:N)=A*bb;
end

figure(3)
mesh(xx,tt,u2)


xlabel('x');
ylabel('t');
title('Solution of Heat equation : Dufort-Frankel')

Error_Dufort =norm(u2-U,Inf)