diff --git a/Exercise_03/NSSC_3.m b/Exercise_03/NSSC_3.m
new file mode 100644
index 0000000..c808100
--- /dev/null
+++ b/Exercise_03/NSSC_3.m
@@ -0,0 +1,43 @@
+%% exercise number 2
+%problem characteristics
+U=0.7;%convection velocity
+delta_x=0.01; %space discretization 
+L=1; %max length of our domain
+T=1; %max time considered
+
+N=L/delta_x; %number of intervals in space
+K=100; %number of intervals in time
+
+%two initial conditions
+square_pulse=@(x) heaviside(x-0.1) - heaviside(x- 0.3); 
+gauss_signal=@(x) exp(-10*(4*x - 1).^2);
+
+[x,t,c] = UW_scheme(L,N,T,K,U,square_pulse);
+ex_sol=@(x,t) square_pulse(x-U*t);
+f1=figure(1);
+for ii=1:K+1
+    clf(f1)
+    hold on 
+    plot(x,c(:,ii)','-bo');
+    plot(x,ex_sol(x,t(ii)));
+    hold off
+    legend("aproximated solution","exact solution")
+    xlim([0 L])
+    ylim([0 1.1])  
+    pause(0.02);
+end
+
+[x,t,c] = UW_scheme(L,N,T,K,U,gauss_signal);
+ex_sol=@(x,t) gauss_signal(x-U*t);
+f2=figure(2);
+for ii=1:K+1
+    clf(f2)
+    hold on 
+    plot(x,c(:,ii)','-bo');
+    plot(x,ex_sol(x,t(ii)));
+    legend("aproximated solution","exact solution")
+    xlim([0 L])
+    ylim([0 1.1])  
+    pause(0.02);
+end
+
diff --git a/Exercise_03/UW_scheme.m b/Exercise_03/UW_scheme.m
new file mode 100644
index 0000000..387dadc
--- /dev/null
+++ b/Exercise_03/UW_scheme.m
@@ -0,0 +1,51 @@
+function [x,t,c] = UW_scheme(xf,N,T,K,U,c0)
+% ----  Numerical solution of the linear advection equation in a periodic domain  ----
+%   c_t+ U * c_x = 0 with domain [0,xf] 
+%   given the initial condition c0.
+% -----------------------------------------------
+% Sintax:
+%   [x,t,c] = UW_scheme(xf,N,T,K,U,c0)
+%
+% Input:
+%   xf end of our domain
+%	N number of space intervals
+%   T max time
+%	K number of time intervals
+%   U convection velocity  
+%   c0 initial condition 
+%
+% Output:
+%   x vector of spatial nodes
+%   t vector of time nodes
+%   c numerical solution 
+
+% Space and time intervals size
+dx=xf/N;
+dt=T/K;
+% initialization of x and t vectors (nodes)
+x=linspace(0,xf,N+1)';
+t=linspace(0,T,K+1)';
+
+% Solution matrix
+c=zeros(N+1,K+1);
+
+% Initial conditions 
+c(:,1) = c0(x);
+
+
+% Creating our matrix
+e = ones(N+1,1);
+B = spdiags([-e,e],[-1,0],N+1,N+1);
+I = speye(N+1);
+
+%"printing" courant number and numerical viscosity
+C0=(U*dt/dx)
+numerical_viscosity= U*dx*(1-C0)/2
+
+%final metrix to compute the solution
+A = I - C0*B;
+
+%finding the solution
+for k=1:K
+    c(1:end,k+1) = A*c(1:end,k);
+end
\ No newline at end of file
diff --git a/Exercise_03/report.tex b/Exercise_03/report.tex
new file mode 100644
index 0000000..cf40f13
--- /dev/null
+++ b/Exercise_03/report.tex
@@ -0,0 +1,298 @@
+\documentclass[a4paper, 10pt]{article}
+
+\usepackage[utf8]{inputenc}
+\usepackage[T1]{fontenc}
+\usepackage{fullpage}
+\usepackage{amsmath, amssymb, amstext, amsthm}
+\usepackage[pdftex]{hyperref}
+\usepackage{xcolor, graphicx}
+\usepackage{tikz}
+\usepackage{listings}
+
+% Document meta into
+\title{NSSC 2 - Assignement 3}
+\author{Bianchi Riccardo, Kapla Daniel, Kuen Jakob, Müller David}
+\date{November 24, 2021}
+% Set PDF title, author and creator.
+\AtBeginDocument{
+    \hypersetup{
+        pdftitle = {NSSC 2 - Exercise 3},
+        pdfauthor = {Group 1},
+        pdfcreator = {\pdftexbanner}
+    }
+}
+
+
+\makeindex
+
+% Setup environments
+% Theorem, Lemma
+\theoremstyle{plain}
+\newtheorem{theorem}{Theorem}
+\newtheorem{lemma}{Lemma}
+\newtheorem{example}{Example}
+% Definition
+\theoremstyle{definition}
+\newtheorem{defn}{Definition}
+% Remark
+\theoremstyle{remark}
+\newtheorem{remark}{Remark}
+
+\DeclareMathOperator*{\argmin}{{arg\,min}}
+\DeclareMathOperator*{\argmax}{{arg\,max}}
+\renewcommand{\t}[1]{{#1^T}}
+\newcommand{\todo}[1]{{\color{red}TODO: #1}}
+
+
+% Default fixed font does not support bold face
+\DeclareFixedFont{\ttb}{T1}{txtt}{bx}{n}{10} % for bold
+\DeclareFixedFont{\ttm}{T1}{txtt}{m}{n}{10}  % for normal
+
+% % Custom colors
+\definecolor{deepblue}{rgb}{0,0,0.5}
+\definecolor{deepred}{rgb}{0.6,0,0}
+\definecolor{deepgreen}{rgb}{0,0.5,0}
+
+% Python style for highlighting
+\newcommand\pythonstyle{\lstset{
+    language=Python,
+    basicstyle=\ttm,
+    morekeywords={self},              % additional keywords
+    keywordstyle=\ttb\color{deepblue},
+    stringstyle=\color{deepgreen},
+    frame=,                           % t, b, or tb for top, bottom env. lines
+    showstringspaces=false,
+    numbers=left,
+    stepnumber=1,
+    numbersep=6pt,
+    numberstyle=\color{gray}\tiny
+}}
+
+% Python environment
+\lstnewenvironment{python}[1][] {
+    \pythonstyle
+    \lstset{#1}
+}{}
+
+\begin{document}
+
+\maketitle
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%                                Exercise 1                                %%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{Numerical solution of the diffusion equation in a finite domain}
+We are given the 1D unsteady diffusion equation
+\begin{equation}\label{eq:ex1}
+    \frac{\partial C}{\partial t} - D \frac{\partial^2 C}{\partial x^2} = 0
+\end{equation}
+with a scalar diffusion coefficient $D = 10^{-6}$. The domain size is $h$ which is discretized with $N_x$ points to give the grid space $\Delta x$ with a time step $\Delta t$. The initial condition is $C = 0$ inside the domain.
+
+Furthermore, denote with $N_t\in\mathbb{N}$ the number of time steps of the simulation where $t_0 = 0$ is the initial time. All $N_t$ discretized time points $t_n = n \Delta t$ leading to the final time point of the simulation as $t_{N_t} = N_t \Delta t$. The discretized $x$ points of the domain are $x_i = (i - 1)\Delta x = (i - 1)h / (N_x - 1)$ for $i = 1, ..., N_x$. This gives exactly $N_x$ grid points with equal distance between adjacent points on the domain $[0, h]$ with the first point $x_1 = 0$ and the last point $x_{N_x} = h$.
+
+In the following we will also use the short hand notation
+\begin{displaymath}
+    C_i^n = C(x_i, t_n), \qquad
+    \partial_x^k C_i^n = \left.\frac{\partial^k C(x, t)}{\partial x^k}\right|_{x = x_i, t = t_n}, \qquad
+    \partial_t^k C_i^n = \left.\frac{\partial^k C(x, t)}{\partial t^k}\right|_{x = x_i, t = t_n}
+\end{displaymath}
+for $k\in\mathbb{N}$ as well as $i$ and $n$ are the space and time discretization indices, respectively.
+
+\subsection{Explicit scheme with Dirichlet/Neumann BC}\label{sec:task01_1}
+Given the Dirichlet boundary conditions $C(0, t) = 0$ and the Neumann BC $\partial_x C(h, t) = 0$. First we derive a $2^{nd}$ order \emph{explicit} finite difference approach with a $1^{st}$ order discretization in time (see Figure~\ref{fig:ex1}).
+
+\begin{figure}[h!]
+    \centering
+    \begin{tikzpicture}[>=latex]
+        \begin{scope}
+            \draw[->] (0.8, 0) -- (5.5, 0) node[anchor = west] {$x$};
+            \draw[->] (1, -0.2) -- (1, 3.5) node[anchor = south] {$t$};
+            \draw[dashed] (0.8, 2) node[anchor = east] {$t_n$} -- (5.5, 2);
+            \draw[dashed] (3, -0.2) node[anchor = north] {$x_i$}
+                -- (3, 3.5) node[anchor = south] {explicit};
+            \foreach \x in {1, ..., 5} {
+                \foreach \y in {0, ..., 3} {
+                    \node[circle, draw, fill = white,
+                        inner sep = 0pt, outer sep = 0pt,
+                        minimum size = 4pt] at (\x, \y) {};
+                }
+            }
+            \node[circle, fill = black,
+                inner sep = 0pt, outer sep = 0pt,
+                minimum size = 4pt] at (3, 2) {};
+            \foreach \x/\y in {2/2, 4/2, 3/3} {
+                \node[circle, fill = gray,
+                    inner sep = 0pt, outer sep = 0pt,
+                    minimum size = 4pt] at (\x, \y) {};
+            }
+        \end{scope}
+        \begin{scope}[xshift = 7cm]
+            \draw[->] (0.8, 0) -- (5.5, 0) node[anchor = west] {$x$};
+            \draw[->] (1, -0.2) -- (1, 3.5) node[anchor = south] {$t$};
+            \draw[dashed] (0.8, 2) node[anchor = east] {$t_n$} -- (5.5, 2);
+            \draw[dashed] (3, -0.2) node[anchor = north] {$x_i$}
+                -- (3, 3.5) node[anchor = south] {implicit};
+            \foreach \x in {1, ..., 5} {
+                \foreach \y in {0, ..., 3} {
+                    \node[circle, draw, fill = white,
+                        inner sep = 0pt, outer sep = 0pt,
+                        minimum size = 4pt] at (\x, \y) {};
+                }
+            }
+            \node[circle, fill = black,
+                inner sep = 0pt, outer sep = 0pt,
+                minimum size = 4pt] at (3, 2) {};
+            \foreach \x/\y in {2/3, 4/3, 3/3} {
+                \node[circle, fill = gray,
+                    inner sep = 0pt, outer sep = 0pt,
+                    minimum size = 4pt] at (\x, \y) {};
+            }
+        \end{scope}
+    \end{tikzpicture}
+    \caption{\label{fig:ex1}Dependency relation for the \emph{explicit} (left) and \emph{implicit} (right) finite difference approximation at a grid point indexed $(i, n)$.}
+\end{figure}
+
+To derive a second order scheme (assuming $C$ is three time continuously differentiable as a function from $[0, h]\times\mathbb{R}^+\to\mathbb{R}$) in space we first considure the Taylor expantion of $C$ with respect to $x$ given by
+\begin{displaymath}
+    C(x + \Delta x, t) = C(x, t) + \partial_x C(x, t)\Delta x + \partial_x^2 C(x, t)\frac{\Delta x^2}{2} + \partial_x^3 C(x, t)\frac{\Delta x^3}{6} + \mathcal{O}(\Delta x^4).
+\end{displaymath}
+Replace $\Delta x$ with $\pm\Delta x$ and add the two equation together, then
+\begin{displaymath}
+    C(x + \Delta x, t) + C(x - \Delta x, t) = 2 C(x, t) + \partial_x^2 C(x, t)\Delta x^2 + \mathcal{O}(\Delta x^4).
+\end{displaymath}
+The first and third order terms drop out due to different signs. Finally, this results in
+\begin{displaymath}
+    \partial_x^2 C(x, t) = \frac{C(x + \Delta x, t) - 2 C(x, t) + C(x - \Delta x, t)}{\Delta x^2} + \mathcal{O}(\Delta x^2)
+\end{displaymath}
+which is an approximation of the second derivative with error proportional to $\Delta x^2$. For the first order approximation in time the same trick can be used except that only the first two terms of the Taylor expantion need to be considured.
+
+Therefore, the $2^{nd}$ order scheme for space and $1^{st}$ order in time at $(i, n)$ using the explicit scheme derives as
+\begin{displaymath}
+    \partial_x^2 \widehat{C}_i^n = \frac{\widehat{C}_{i-1}^{n} - 2\widehat{C}_{i}^{n} + \widehat{C}_{i+1}^{n}}{\Delta x^2},
+    \qquad
+    \partial_t \widehat{C}_i^n = \frac{\widehat{C}_{i}^{n+1} - \widehat{C}_{i}^{n}}{\Delta t}.
+\end{displaymath}
+Substitution into \eqref{eq:ex1} yields after rearranging the update rule for internal points as
+\begin{equation}\label{eq:task01_1_update}
+    \widehat{C}_i^{n+1} = \widehat{C}_i^n + \frac{D\Delta t}{\Delta x^2}(\widehat{C}_{i-1}^{n} - 2\widehat{C}_{i}^{n} + \widehat{C}_{i+1}^{n}).
+\end{equation}
+This holds at $x_i$ where $i = 2, ..., N_x - 1$, or in other words everywhere inside the space domain excluding the boundary. The boundary needs to be handled seperately. The left boundary condition is a Dirichlet constraint which is simply constanct giving $C_1^n = 0$ for all time while the Neumann condition on the right requires an additional discretization step for computing the next value. Therefore, we take the second order Taylor expantion of $C$ with respect to $x$ in the negative direction given by
+\begin{displaymath}
+    C(x - \Delta x, t) = C(x, t) - \partial_x C(x, t)\Delta x + \mathcal{O}(\Delta x^2).
+\end{displaymath}
+As we are interesetd in the boundary value for the Neumann BC, the derivative at $x = h$ is given which leads to
+\begin{displaymath}
+    C(h, t) = C(h - \Delta x, t) + \partial_x C(h, t)\Delta x + \mathcal{O}(\Delta x^2)
+\end{displaymath}
+which is second order accurate. Therefore, the right boundary value for our boundary condition $\partial_x C(h, t_n) = 0$ at time $t_n$ is
+\begin{displaymath}
+    \widehat{C}_{N_x}^n = \widehat{C}_{N_x - 1}^n.
+\end{displaymath}
+See \texttt{task01\_1-2.py} or \texttt{Mixed\_EA.m} for an implementation of this scheme.
+
+\subsubsection{Stability} \todo{Prove that the update scheme is unstable for $d = \frac{D\Delta x^2}{\Delta t} > 0.5$}
+
+\subsection{Explicit scheme with Dirichlet BC at both boundaries}\label{sec:task01_2}
+Not we have the same setting as in Section~\ref{sec:task01_1} except for Dirichlet boundary conditions on both sides given by $C(0, t) = 1$ and $C(h, t) = 0$. This means tha the left and right values are constants and internal nodes follow the same update scheme as in \eqref{eq:task01_1_update}.
+
+\todo{include pictures and compare with Section~\ref{sec:task01_1}}
+
+\subsection{Implicit scheme with Dirichlet/Neumann BC}\label{sec:task01_3}
+The following computations are similar in nature to the explicit scheme, therefore well keep it short. The implicit discretization for the derivatives is (see: Figure~\ref{fig:ex1}) then
+\begin{displaymath}
+    \partial_x^2 C_{i}^{n+1} = \frac{\widehat{C}_{i-1}^{n+1} - 2\widehat{C}_{i}^{n+1} + \widehat{C}_{i+1}^{n+1}}{\Delta x^2},
+    \qquad
+    \partial_t \widehat{C}_i^n = \frac{\widehat{C}_{i}^{n+1} - \widehat{C}_{i}^{n}}{\Delta t}.
+\end{displaymath}
+Substitution into \eqref{eq:ex1} gives the implicit (inverse) update rule
+\begin{displaymath}
+    \widehat{C}_i^{n} = -d \widehat{C}_{i-1}^{n+1} + (1 + 2 d)  \widehat{C}_i^{n+1} - d \widehat{C}_{i+1}^{n+1}
+\end{displaymath}
+for $i = 2, ..., N_x - 1$ and $d = \frac{D \Delta x^2}{\Delta t}$. The boundary conditions are ether $\widehat{C}_1^{n+1} = 0$ for the left Dirichlet and $\widehat{C}_{N_x}^{n+1} = \widehat{C}_{N_x-1}^{n+1}$ as the right Neumann BC. By collecting all coefficients of the $i = 2, ..., N_x - 1$ into a $(N_x - 2)\times (N_x - 2)$ trigiagonal matrix
+\begin{displaymath}
+    A = \begin{pmatrix}
+        1+2d & -d &   \\
+        -d & 1+2d & -d &   \\
+          & -d & 1+2d & \ddots \\
+          &   & \ddots & \ddots & -d \\
+          & & & -d & 1+2d & -d \\
+          & & &   & -d & 1+d
+      \end{pmatrix}
+\end{displaymath}
+\todo{fix $1+d$ which is first order accurate!!!}
+
+where the last entry $A_{N_x, N_x} = 1 - d$ is due to the Neumann BC $\widehat{C}_{N_x}^{n+1} = \widehat{C}_{N_x-1}^{n+1}$ which means
+\begin{displaymath}
+    \widehat{C}_i^{N_x-1} = -d \widehat{C}_i^{N_x-2} + (1 + 2 d) \widehat{C}_i^{N_x - 1} - d \widehat{C}_i^{N_x}
+        = -d \widehat{C}_i^{N_x-2} + (1 + d) \widehat{C}_i^{N_x - 1}.
+\end{displaymath}
+Finaly, we end up with the implicit update rule
+\begin{displaymath}
+    \widehat{C}^{n+1} = A \widehat{C}^{n}
+\end{displaymath}
+which is perfomed by solving the linear system for $\widehat{C}^{n+1}$.
+
+\todo{include pictures and compare with Section~\ref{sec:task01_1}}
+
+\subsection{Second order in time for implicit scheme with Dirichlet/Neumann BC}
+To derive the second order accurate scheme in time we employ the Crank-Nicolson approach to derive the descritization at the time midpoints $t_{n+1/2} = (n + 1/2)\Delta t$ which we index with $n+1/2$.
+
+The second order im space is identical except the evaluation points
+\begin{displaymath}
+    \partial_x^2 C_i^{n+1/2} = \frac{C_{i-1}^{n+1/2} - 2C_{i}^{n+1/2} + C_{i+1}^{n+1/2}}{\Delta x^2} + \mathcal{O}(\Delta x^2)
+\end{displaymath}
+and for the time descritization we take the Taylor expantion
+
+
+\begin{displaymath}
+    \partial_{t}C_i^{n+1/2} = \frac{C_i^{n+1} - C_i^{n}}{2\Delta t} + \mathcal{O}(\Delta t^2)
+\end{displaymath}
+\todo{check if this follows from $\exists \partial_x^4 C$ beeing continuous}
+\begin{displaymath}
+    \partial_x^2 C_i^{n+1/2} = \frac{\partial_x^2 C_i^{n} + \partial_x^2 C_i^{n+1}}{2} + \mathcal{O}(\Delta x^2)
+\end{displaymath}
+
+Meaning
+\begin{align*}
+    0 = \partial_{t}C_i^{n+1/2}-D\partial_x^2 C_i^{n+1/2} &= \frac{C_i^{n+1} - C_i^{n}}{2\Delta t}
+     -D\frac{\partial_x^2 C_i^{n} + \partial_x^2 C_i^{n+1}}{2} + \mathcal{O}(\Delta t^2 + ?) \\
+     &= \frac{C_i^{n+1} - C_i^{n}}{2\Delta t} -
+         D\frac{C_{i-1}^{n} + C_{i-1}^{n+1} - 2(C_{i}^{n} + C_{i}^{n+1}) + C_{i+1}^{n} + C_{i+1}^{n+1}}{2\Delta x^2}
+\end{align*}
+rearranging yields
+\begin{displaymath}
+    -d C_{i-1}^{n+1} + (1 + 2d) C_i^{n+1}-dC_{i+1}^{n+1}
+        = d C_{i-1}^n + (1 - 2d)C_i^n + dC_{i+1}^n
+        = -(-d C_{i-1}^n + (1 + 2d)C_i^n - dC_{i+1}^n) + 2C_i^n
+\end{displaymath}
+Now with second order approximation of the Neumann BC
+\begin{displaymath}
+    0\overset{!}{=} \partial_x C_{N_x}^{n+1} = \frac{C_{N_x-2}^{n+1} - C_{N_x}^{n+1}}{2\Delta x} + \mathcal{O}(\Delta x^2)
+\end{displaymath}
+which means that (\todo{check why?! should be at $N_x - 1$}) for index $i = N_x - 1$ 
+\begin{align*}
+    -d C_{N_x-2}^{n+1} + (1 + 2d) C_{N_x-1}^{n+1}-dC_{N_x}^{n+1}
+        &= -(-d C_{N_x-2}^n + (1 + 2d)C_{N_x-1}^n - dC_{N_x}^n) + 2C_{N_x-1}^n \\
+    -d C_{N_x-2}^{n+1} + (1 + 2d) C_{N_x-1}^{n+1}-d C_{N_x-2}^{n+1}
+        &= -(-d C_{N_x-2}^n + (1 + 2d)C_{N_x-1}^n - d C_{N_x-2}^{n}) + 2C_{N_x-1}^n \\
+    -2d C_{N_x-2}^{n+1} + (1 + 2d) C_{N_x-1}^{n+1}
+        &= -(-2d C_{N_x-2}^n + (1 + 2d)C_{N_x-1}^n) + 2C_{N_x-1}^n
+\end{align*}
+
+\begin{displaymath}
+    A C^{n+1} = (2 I - A) C^n
+\end{displaymath}
+with 
+\begin{displaymath}
+    A = \begin{pmatrix}
+        1 & 0 & 0 \\
+        -d & 1+2d & -d \\
+        & -d & 1+2d & -d \\
+        & & -d & 1+2d & \ddots \\
+        & & & \ddots & \ddots & -d \\
+        & & & & -d & 1+2d & -d \\
+        & & & & & -2d & 1+2d
+    \end{pmatrix}
+\end{displaymath}
+
+\end{document}
diff --git a/Exercise_03/task01_3.py b/Exercise_03/task01_3.py
index b824a47..021624d 100644
--- a/Exercise_03/task01_3.py
+++ b/Exercise_03/task01_3.py
@@ -39,7 +39,7 @@ for t in range(nt):
         plt.ylim([0, 1.2])
         plt.savefig(f"plots/task01_3_{i:0>5}.png")
         i += 1
-    # update solution using the explicit schema
+    # update solution using the implicit schema
     C = np.linalg.solve(T, C)
     # fix BC conditions (theoretically, they are set by the update but for
     # stability reasons (numerical) we enforce the correct values)
diff --git a/Exercise_03/task01_4.py b/Exercise_03/task01_4.py
new file mode 100644
index 0000000..68927b5
--- /dev/null
+++ b/Exercise_03/task01_4.py
@@ -0,0 +1,74 @@
+# Task 1.4
+import numpy as np
+import scipy, scipy.sparse, scipy.linalg
+from matplotlib import pyplot as plt
+
+# plotting helper function
+def plot(C, x, t, name):
+    plt.clf()
+    plt.plot(x, C)
+    plt.xlim([0, x[-1]])
+    plt.ylim([0, 1.2])
+    plt.title(f"{name} - time: {t: >12.0f}")
+    plt.savefig(f"plots/{name}_{plot.fignr:0>5}.png")
+    plot.fignr += 1
+
+# set static variable as function attribute
+plot.fignr = 1
+# and initialize a figure
+plt.figure(figsize = (8, 6), dpi = 100)
+
+# Config
+D = 1e-6            # diffusion coefficient
+h = 1               # space domain (max x size)
+T = 1e5             # solution end time
+nx = 50             # nr of space discretization points
+nt = 2000           # nr of time points
+
+# derived constants
+dx = h / (nx - 1)   # space step size
+dt = T / (nt - 1)   # time step size
+d = dt * D / dx**2  # stability/stepsize coefficient
+
+# setup matrices for implicit update scheme `A C^n+1 = (2 I - A) C^n = B C^n`
+# The matrix `A` is represented as a `3 x (nx - 1)` matrix as expected by
+# `scipy.linalg.solve_banded` representing the 3 diagonals of the tridiagonal
+# matrix `A`.
+A = np.zeros((3, nx - 1))
+# upper minor diagonal
+A[0, 2:] = -d
+# main diagonal
+A[1, 0] = 1
+A[1, 1:] = 1 + 2 * d
+# and lower minor diagonal
+A[2, :-2] = -d
+A[2, -2] = -2 * d
+# The right hand side matrix `B = 2 I - A` which is represented as a sparse
+# matrix in diag. layout which allows for fast matrix vector multiplication
+B = -A
+B[1, ] += 2
+B = scipy.sparse.spdiags(B, (1, 0, -1), B.shape[1], B.shape[1])
+
+# set descritized space evaluation points `x` of `C`
+x = np.linspace(0, h, nx)
+
+# Set initial solution
+C = np.zeros(nx)
+C[0] = 1
+C[-1] = C[-3] # (0 = 0)
+
+for n in range(nt):
+    # generate roughly 400 plots
+    if n % (nt // 400) == 0:
+        plot(C, x, n * dt, "task01_4")
+    # update solution using the implicit schema
+    C[:-1] = scipy.linalg.solve_banded((1, 1), A, B @ C[:-1])
+    # set/enforce boundary conditions
+    C[0] = 1            # left Dirichlet (theoretically not needed)
+    C[-1] = C[-3]       # right Neumann condition
+
+# plot final state
+plot(C, x, T, "task01_4")
+
+# to convert generated image sequence to video use:
+# $> ffmpeg -r 60 -i plots/task01_4_%05d.png -pix_fmt yuv420p video_1_4.mp4
diff --git a/Exercise_03/video_1_4.mp4 b/Exercise_03/video_1_4.mp4
new file mode 100644
index 0000000..ab99238
Binary files /dev/null and b/Exercise_03/video_1_4.mp4 differ