diff --git a/Exercise_03/Dirichlet_EA.m b/Exercise_03/Dirichlet_EA.m
new file mode 100755
index 0000000..1243748
--- /dev/null
+++ b/Exercise_03/Dirichlet_EA.m
@@ -0,0 +1,59 @@
+function [x,t,u] = Dirichlet_EA(L,N,T,K,c1,c2,f,u0,D)
+% ----  Solution of heat equation ----
+%   u_t - u_xx = f in the interval [-L,L] (doesn't matter if we change it)
+%   Dirichlet BC
+%  
+% -----------------------------------------------
+% Sintassi:
+%   [x,t,u]=calore_template(L,N,T,K,c1,c2,fun,u0)
+%
+% Input:
+%   L Half of the width (-L,L) 
+%	N number of intervals in (-L,L)
+%   T max time (0,T)
+%	K number of intervals in (0,T)
+%   c1 Dirichlet BC in x=-L
+%   c2 Dirichlet BC in x=L
+%   f force function
+%   u0 initial condition in t=0
+%
+% Output:
+%   x vector of the spatial nodes
+%   t vector of the time nodes
+%   u numeric solution of the problem
+
+% discretisation step in time and space
+h=L/N; %space
+tau=T/K; %time
+% Initialization of t
+t=linspace(0,T,K+1)';
+% initialization of x
+x=linspace(0,L,N+1);
+
+% Initialization of the matrix solution u
+u=zeros(N+1,K+1);
+% Intial conditions
+u(:,1)=u0(x);
+
+% BC
+u(1,:)=c1(t);
+u(end,:)=c2(t);
+
+% Creation of the matrix A
+% We obtained this matrix using the Central Discretization method
+e=ones(N-1,1);
+A=spdiags([-e,2*e,-e],[-1,0,1],N-1,N-1)/(h^2);
+I=speye(N-1,N-1);
+
+% Compute the solution
+% What we are doing here is to compute the solution in the interval for
+% each time k
+for k=1:K
+    % Assembly of the force
+    F=f(x(2:end-1),t(k));
+    % Correction of the force using the BC
+    F(1)=F(1) + D*c1(t(k))/(h^2);
+    F(end)=F(end) + D*c2(t(k))/(h^2);
+    % Solution using Eulero forward
+    u(2:end-1,k+1) = ((I - tau*D*A)*u(2:end-1,k))' + tau*F;
+end
diff --git a/Exercise_03/Mixed_EA.m b/Exercise_03/Mixed_EA.m
new file mode 100755
index 0000000..1ac4aee
--- /dev/null
+++ b/Exercise_03/Mixed_EA.m
@@ -0,0 +1,63 @@
+function [x,t,u] = Mixed_EA(L,N,T,K,c1,c2,f,u0,D)
+% ----  Solution of heat equation ----
+%   u_t - u_xx = f in the interval [-L,L] (doesn't matter if we change it)
+%   Dirichlet BC
+%  
+% -----------------------------------------------
+% Sintassi:
+%   [x,t,u]=calore_template(L,N,T,K,c1,c2,fun,u0)
+%
+% Input:
+%   L Half of the width (-L,L) 
+%	N number of intervals in (-L,L)
+%   T max time (0,T)
+%	K number of intervals in (0,T)
+%   c1 Dirichlet BC in x=-L
+%   c2 Dirichlet BC in x=L
+%   f force function
+%   u0 initial condition in t=0
+%
+% Output:
+%   x vector of the spatial nodes
+%   t vector of the time nodes
+%   u numeric solution of the problem
+
+% discretisation step in time and space
+h=L/N; %space
+tau=T/K; %time
+% Initialization of t
+t=linspace(0,T,K+1)';
+% initialization of x
+x=linspace(0,L,N+1);
+
+% Initialization of the matrix solution u
+u=zeros(N+1,K+1);
+% Intial conditions
+u(:,1)=u0(x);
+
+% BC Dirichlet
+u(1,:)=c1(t);
+
+% Creation of the matrix A
+% We obtained this matrix using the Central Discretization method
+e=ones(N-1,1);
+A=spdiags([-e,2*e,-e],[-1,0,1],N-1,N-1)/(h^2);
+%modify the matrix such that we can use the noimann condition to calculate
+%the last node.
+%In this case we want to use a second order decentralized approximation 
+A(end,end-1)=-2/(h^2);
+I=speye(N-1,N-1);
+
+% Compute the solution
+% What we are doing here is to compute the solution in the interval for
+% each time k
+for k=1:K
+    % Assembly of the force
+    F=f(x(2:end-1),t(k));
+    % Correction of the force using the BC
+    F(1)=F(1) + D*c1(t(k))/(h^2);
+    F(end)=F(end) + 2*h*c2;
+    % Solution using Eulero forward
+    u(2:end-1,k+1) = ((I - D*tau*A)*u(2:end-1,k))' + tau*F;
+    u(end,k+1)=u(end-2,k+1) + c2*2*h*D;
+end
diff --git a/Exercise_03/Mixed_EI.m b/Exercise_03/Mixed_EI.m
new file mode 100755
index 0000000..635f458
--- /dev/null
+++ b/Exercise_03/Mixed_EI.m
@@ -0,0 +1,62 @@
+function [x,t,u] = Mixed_EI(L,N,T,K,c1,c2,f,u0,D)
+% ----  Risoluzione dell'equazione del calore ----
+%   u_t - u_xx = f nell'intervallo [-L,L] 
+%   con condizioni al bordo di Dirichlet
+%   e condizioni iniziali.
+% -----------------------------------------------
+% Sintassi:
+%   [x,t,u]=calore_template(L,N,T,K,c1,c2,fun,u0)
+%
+% Input:
+%   L semiampiezza intervallo spaziale (-L,L) 
+%	N numero di sottointervalli in (-L,L)
+%   T estremo finale intervallo temporale (0,T)
+%	K numero di sottointervalli in (0,T)
+%   c1 funzione che descrive la condizione di Dirichlet in x=-L
+%   c2 funzione che descrive la condizione di Dirichlet in x=L
+%   f funzione che descrive il termine noto dell'equazione
+%   u0 funzione che descrive la condizione iniziale in t=0
+%
+% Output:
+%   x vettore dei nodi spaziali
+%   t vettore dei nodi temporali
+%   u soluzione numerica 
+%     del problema
+
+% Calcolo passo di discretizzazione in spazio e tempo
+h=L/N;
+tau=T/K;
+% Inizializzazione del vettore t
+t=linspace(0,T,K+1);
+% Inizializzazione del vettore x
+x=linspace(0,L,N+1);
+
+% Inizializzazione della matrice soluzione u
+u=zeros(N+1,K+1);
+% Condizione iniziale
+u(:,1)=u0(x);
+
+% BC Dirichlet
+u(1,:)=c1(t);
+
+% Creation of the matrix A
+% We obtained this matrix using the Central Discretization method
+e=ones(N-1,1);
+A=spdiags([-e,2*e,-e],[-1,0,1],N-1,N-1)/(h^2);
+%modify the matrix such that we can use the noimann condition to calculate
+%the last node.
+%In this case we want to use a second order decentralized approximation 
+A(end,end-1)=-2/(h^2);
+I=speye(N-1,N-1);
+
+
+for k=1:K
+    % Assemblaggio termine noto
+    F=f(x(2:end-1),t(k+1));
+    % Correzione del termine noto con le condizioni al bordo
+    F(1)=F(1) + c1(t(k+1))/(h^2);
+    F(end)=c2*h*2;
+    % Risoluzione del problema
+    u(2:end-1,k+1) = ((I+tau*A)\u(2:end-1,k))' + tau*F;
+    u(end,k+1)=u(end-2,k+1) + c2*2*h;
+end
diff --git a/Exercise_03/NSSC_1.m b/Exercise_03/NSSC_1.m
new file mode 100755
index 0000000..f548749
--- /dev/null
+++ b/Exercise_03/NSSC_1.m
@@ -0,0 +1,93 @@
+%% eulero forward DIrichlet
+L=pi;
+T=2;
+f=@(x,t) 0*x.*t;
+c1=@(t) 1+0*t;
+c2=@(t) 0*t;
+u0=@(x) 0*x;
+D=0.5;
+%uex=@(x,t) cos(x).*exp(t);
+
+N=10;
+K=200;
+[x,t,u]=Dirichlet_EA(L,N,T,K,c1,c2,f,u0,D);
+
+figure(1)
+for ii=1:K+1 
+    plot(x,u(:,ii)');
+    xlim([0 L])
+    pause(0.05);
+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 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)