% 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 explicit Euler 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

% set up grid in x and t
N=20;
x=linspace(0,1,N+1);
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;

% explicit Euler loop
for m=1:M
    U(1)=2*pi*exp(-4*pi^2*t(m+1));
    for j=2:N
        U(j)=Uold(j)+mu*(Uold(j+1)-2*Uold(j)+Uold(j-1));
    end
    U(N+1)=2*pi*exp(-4*pi^2*t(m+1));
    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)))