% Solve 1D heat equation on [0,1]x[0,T] with initial condition
% u(x,0)=sin(pi*x)+2*pi*cos(2*pi*x) and boundary conditions
% u(0,t)=2*pi*exp(-4*pi^2*t)=u(1,t). The exact solution is
% u(x,t)=exp(-pi^2*t)*sin(pi*x)+2*pi*exp(-4*pi^2*t)*cos(2*pi*x).
% We use central differences in space and the theta method for timestepping.

% exact solution
uexact=@(x,t) exp(-pi^2*t)*sin(pi*x)+2*pi*exp(-4*pi^2*t)*cos(2*pi*x);
figure(1), clf
fplot(@(x) uexact(x,0.1),[0,1],'b','LineWidth',2), hold on

theta=1;

% set up grid in x and t
N=40;
x=linspace(0,1,N+1); x=x(:);
dx=x(2)-x(1);
dt=dx^2/4;
t=0:dt:0.1;
mu=dt/dx^2;
M=length(t)-1;

% initial condition
U=sin(pi*x)+2*pi*cos(2*pi*x);
Uold=U;

% theta method timestepping
% first set up matrices which are independent of time
e=ones(N+1,1);
A=spdiags([e -2*e e], -1:1, N+1,N+1);
A(1,:)=0; A(N+1,:)=0;
LHSmat=speye(N+1)-mu*theta*A;
Ip=speye(N+1); Ip(1,1)=0; Ip(N+1,N+1)=0;
RHSmat=Ip+mu*(1-theta)*A;
g=zeros(N+1,1);

% now loop over time
for m=1:M
    g(1)=2*pi*exp(-4*pi^2*t(m+1));
    g(N+1)=g(1);
    U=LHSmat\(RHSmat*Uold+g);
    Uold=U;
end

plot(x,U,'r','LineWidth',2)
xlabel('x'), ylabel('u(x,0.1)'), legend('Exact solution','Numerical solution')

figure(2)
plot(x,U-uexact(x,0.1))
max(abs(U-uexact(x,0.1)))