diff --git a/.DS_Store b/.DS_Store
new file mode 100644
index 0000000..dbf098b
Binary files /dev/null and b/.DS_Store differ
diff --git a/Exercise_03/NSSC_1.m b/Exercise_03/NSSC_1.m
new file mode 100755
index 0000000..e694ce7
--- /dev/null
+++ b/Exercise_03/NSSC_1.m
@@ -0,0 +1,107 @@
+%% Euler forward Dirichlet
+L = 1;  % domain size (in the lecture notes this is denote h)
+T = 2;  % time limit (max time)
+f = @(x,t) 0*x.*t;    % rhs of the more general equation `u_t - d u_xx = f`
+c1 = @(t) 1+0*t;      % _right_ boundary condition
+c2 = @(t) 0*t;        % _left_ boundary condition
+u0 = @(x) 0*x;        % initial values
+D = 0.5;              % diffusion parameter `d` in `u_t - d u_xx = f`
+%uex = @(x,t) cos(x).*exp(t);
+
+N = 10;               % nr. of _space_ discretization points
+K = 200;              % nr. of _time_ discretization points
+[x, t, u] = Dirichlet_EA(L, N, T, K, c1, c2, f, u0, D);
+
+% Report stability condition `D Delta T / (Delta x)^2 > 0.5`
+Delta_T = T / K;
+Delta_x = L / N;
+d = D * Delta_T / Delta_x^2;
+fprintf("Stability Condition: 0.5 >= D * Delta_T / Delta_x^2 = %f\n", d)
+if d > 0.5
+    fprintf("-> NOT Stable\n")
+else
+    fprintf("-> Stable\n")
+end
+
+figure(1)
+for ii = 1:K+1 % iterates time
+    hold on
+    plot(x, u(:, ii)');
+    xlim([0 L])
+    pause(0.05);
+    hold off
+end
+
+% 3D plot of space solution over time
+space = linspace(0,L,101);
+time = linspace(0,T,201);
+[xx,yy] = meshgrid(time,space);
+%exsol = uex(yy,xx);
+figure(2)
+mesh(t,x,u)
+%figure(2)
+%mesh(xx,yy,exsol)
+
+%% Eulero forward Mixed BC
+
+L=2*pi;
+T=5;
+f=@(x,t) 0*x.*t;
+c1=@(t) 1+0*t;
+c2=0;
+u0=@(x) 0*x;
+D=1.1;
+%uex=@(x,t) cos(x).*exp(t);
+
+N=25;
+K=200;
+[x,t,u]=Mixed_EA(L,N,T,K,c1,c2,f,u0,D);
+
+figure(1)
+for ii=1:K+1 
+    plot(x,u(:,ii)');
+    xlim([0 L])
+    ylim([0 1.5])
+    pause(0.02);
+end
+
+space=linspace(0,L,101);
+time=linspace(0,T,201);
+[xx,yy]=meshgrid(time,space);
+%exsol=uex(yy,xx);
+figure(2)
+mesh(t,x,u)
+%figure(2)
+%mesh(xx,yy,exsol)
+
+%% Eulero Backward Mixed BC
+
+L=2*pi;
+T=5;
+f=@(x,t) 0*x.*t;
+c1=@(t) 1+0*t;
+c2=0;
+u0=@(x) 0*x;
+D=1;
+%uex=@(x,t) cos(x).*exp(t);
+
+N=25;
+K=200;
+[x,t,u]=Mixed_EI(L,N,T,K,c1,c2,f,u0,D);
+
+figure(1)
+for ii=1:K+1 
+    plot(x,u(:,ii)');
+    xlim([0 L])
+    ylim([0 1.5])
+    pause(0.02);
+end
+
+space=linspace(0,L,101);
+time=linspace(0,T,201);
+[xx,yy]=meshgrid(time,space);
+%exsol=uex(yy,xx);
+figure(2)
+mesh(t,x,u)
+%figure(2)
+%mesh(xx,yy,exsol)
diff --git a/Exercise_04/Zienkiewicz_Chapter_5.pdf b/Exercise_04/Zienkiewicz_Chapter_5.pdf
new file mode 100644
index 0000000..d3e029e
--- /dev/null
+++ b/Exercise_04/Zienkiewicz_Chapter_5.pdf
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:d3430c289a2cf72964e61b1a402d063473b1d866cfbdfe744c1b8a18949ff4a9
+size 974296
diff --git a/Exercise_04/inputfile_group_1.txt b/Exercise_04/inputfile_group_1.txt
new file mode 100644
index 0000000..217e2f7
--- /dev/null
+++ b/Exercise_04/inputfile_group_1.txt
@@ -0,0 +1,39 @@
+"""
+NSSCII - FEM.
+Input to be parsed through.
+
+SI-units to be used:
+    + T in K
+    + L in m
+    + k in W/(mK)
+    + q in W/m^2 - ad Neumann
+    + P in W - ad nodal forces
+"""
+
+# Group number.
+groupnr = 1
+
+# Length in x- and y-direction.
+L = 0.01
+
+# Thickness (z-direction).
+hz = 0.0005
+
+# Thermal conductivity (k=k_xx=k_yy, k_xy = 0.).
+k = 429.
+
+# Factor c for modifying thermal conductivity k for
+# elements in elements_to_be_modified.
+c = 10.
+
+# Elements to be modified.
+elements_to_be_modified = [
+                          41-47,
+                          59-63,
+                          77-79,
+                          95
+                          ]
+
+# Boundary conditions.
+q(y=0) = 2000000.
+T(y=L) = 293.
diff --git a/Exercise_04/nsscii-exc4-fem.pdf b/Exercise_04/nsscii-exc4-fem.pdf
new file mode 100644
index 0000000..f638419
--- /dev/null
+++ b/Exercise_04/nsscii-exc4-fem.pdf
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:0dd52c65911bca256fa5565f888fcfb1485ae49abe0525dc256f4d0c2bb50a83
+size 202797
diff --git a/Exercise_04/print_HTP.py b/Exercise_04/print_HTP.py
new file mode 100644
index 0000000..b969c06
--- /dev/null
+++ b/Exercise_04/print_HTP.py
@@ -0,0 +1,36 @@
+def print_HTP(H, T, P, filename="output.txt"):
+    """
+    Print matrices to .txt-file (name of file = filename).
+    H... overall assembled stiffness matrix
+    T... nodal temperature vector
+    P... nodal force vector
+
+    Make sure, that your system of equations is sorted by
+    ascending node numbers, i.e., N1 N2 ... N100.
+    """
+
+    F = open(filename, 'w')
+
+    F.write("Stiffness matrix H: \n")
+    for row in H:
+        for col in row:
+            outline = "{0:+8.4e},".format(col)
+            F.write("{0:11s}".format(str(outline)))
+        F.write("\n")
+
+    F.write("Temperature T: \n")
+    for row in T:
+        for col in row:
+            outline = "{0:+8.4e},".format(col)
+            F.write("{0:11s} \n".format(str(outline)))
+
+
+    F.write("Force vector P: \n")
+    for row in P:
+        for col in row:
+            outline = "{0:+8.4e},".format(col)
+            F.write("{0:11s} \n".format(str(outline)))
+
+    F.close()
+
+    return None