
clear

f = @(x) 0.1*sin(0.5*pi*x);

h = 0.2;
X = (0:h:2)';
n = length(X);

U0 = max(1-5*abs(1-X),0);
Uee = U0;
Ucn = U0;

M = 1/h^2*(-2*diag(ones(n-2,1))+diag(ones(n-3,1),1)+diag(ones(n-3,1),-1));

plot(X,U0)

k = 0.005;
k = 0.025;
k = 0.1;
T = 0:k:0.5;

A = (eye(9,9)-k/2*M)^(-1);
for i=1:length(T)-1
    Uee(2:end-1,i+1) = (eye(9,9)+k*M)*Uee(2:end-1,i)+k*f(X(2:end-1));
    Ucn(2:end-1,i+1) = A*((eye(9,9)+k/2*M)*Ucn(2:end-1,i)+k*f(X(2:end-1)));
end

figure(1)
mesh(T,X,Uee)
title(['Expl. Euler, $h=',num2str(h),',\;k=',num2str(k),'$'],'interpreter','latex')
axis tight
set(gca,'fontsize',15)

figure(2)
mesh(T,X,Ucn)
title(['Crank-Nicolson, $h=',num2str(h),',\;k=',num2str(k),'$'],'interpreter','latex')
axis tight
set(gca,'fontsize',15)

k = 0.05;
k = 0.1;
T = 0:k:0.5;
Uie = U0;
for i=1:length(T)-1
    Uie(2:end-1,i+1) = (eye(9,9)-k*M)\(Uie(2:end-1,i)+k*f(X(2:end-1)));
end

figure(3)
mesh(T,X,Uie)
title(['Impl. Euler-Verfahren, $h=',num2str(h),',\;k=',num2str(k),'$'],'interpreter','latex')
axis tight
set(gca,'fontsize',15)
