diff --git a/symb_abs.m b/symb_abs.m
new file mode 100644
index 0000000000000000000000000000000000000000..d896d4eec55392b32666b640a103d511c84d8ac4
--- /dev/null
+++ b/symb_abs.m
@@ -0,0 +1,100 @@
+%% Tutorial example on symbolic control
+%% Author: Antoine GIRARD, Université Paris-Saclay, CNRS, CentraleSupélec, Laboratoire des signaux et systèmes, 91190, Gif-sur-Yvette, France.
+%% Date: 2024
+
+% Computation of a symbolic model for the system defined in dynamics.m
+
+clear 
+
+dynamics;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Abstraction
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+tic
+
+% This part of the script is dedicated to computing the symbolic model
+
+% Space discretization
+d_x=(bound_x(:,2)-bound_x(:,1))./n_x;    
+p_x=[1;cumprod(n_x)];        
+n_xi=prod(n_x);                
+grid_xi=reshape(1:n_xi,n_x');
+
+% abstraction function
+q=@(x) coord2state(floor((x-bound_x(:,1))./d_x)+1,p_x);
+
+% Input discretization
+d_u=(bound_u(:,2)-bound_u(:,1))./(n_u-1);
+p_u=[1;cumprod(n_u)];        
+n_sigma=prod(n_u);
+
+% concretisation function
+val_u=zeros(length(bound_u),n_sigma);
+for i=1:n_sigma
+    val_u(:,i)=bound_u(:,1)+(state2coord(i,p_u)-[1;1]).*d_u;
+end
+
+p=@(sigma) val_u(:,sigma);
+
+% Disturbance
+d_w=(bound_w(:,2)-bound_w(:,1));   
+w_c=0.5*(bound_w(:,1)+bound_w(:,2)); 
+
+disp('Computation of the symbolic model')
+h = waitbar(0,'Abstraction in progress...');
+
+% Computation of the transition relation
+% g stores the lower and upper successor states 
+% for each state under the action of each input
+
+g=zeros(n_xi,n_sigma,2);
+robust_margins=d_x;
+for xi=1:n_xi 
+    
+    waitbar(xi/n_xi,h)
+    
+    c_xi=state2coord(xi,p_x);         
+    x_c=bound_x(:,1)+(c_xi-0.5).*d_x; % Coordinates of the center of the cell
+    
+    % Compute the transition relation 
+    for sigma=1:n_sigma
+        
+        % Over approximation of the reachable set
+        x_c_succ=f(x_c, val_u(:,sigma), w_c); % compute the image of the center point
+        d_x_succ=0.5*Jf_x(val_u(:,sigma))*d_x+0.5*Jf_w(val_u(:,sigma))*d_w; % bound the deviation 
+        reach=[x_c_succ-d_x_succ,x_c_succ+d_x_succ]; % reachable set
+        
+        %Special treatment of the third coordinate (angle) 
+        %not needed in general
+        if reach(3,1)<-pi && reach(3,2)>=-pi
+            reach(3,1)=reach(3,1)+2*pi;
+        elseif reach(3,1)<-pi && reach(3,2)<-pi
+            reach(3,1)=reach(3,1)+2*pi;
+            reach(3,2)=reach(3,2)+2*pi; 
+        elseif reach(3,2)>pi && reach(3,1)<=pi 
+            reach(3,2)=reach(3,2)-2*pi;
+        elseif reach(3,2)>pi && reach(3,1)>pi
+            reach(3,2)=reach(3,2)-2*pi;
+            reach(3,1)=reach(3,1)-2*pi;
+        end
+        
+        % Store the transition relation
+        if all(reach(:,1)>=bound_x(:,1)) && all(reach(:,2)<=bound_x(:,2))
+            min_succ=coord2state(floor((reach(:,1)-bound_x(:,1))./d_x)+1,p_x); 
+            max_succ=coord2state(ceil((reach(:,2)-bound_x(:,1))./d_x),p_x);
+            g(xi,sigma,:)=[min_succ,max_succ];
+            robust_margins=min(robust_margins,...
+                           (reach(:,1)-bound_x(:,1))-d_x.*floor((reach(:,1)-bound_x(:,1))./d_x));
+            robust_margins=min(robust_margins,...
+                           (-reach(:,2)+bound_x(:,1))+d_x.*ceil((reach(:,2)-bound_x(:,1))./d_x));           
+        end
+        % g is 0 if and only if reach set outside the state constraints    
+    end
+
+end
+close(h) 
+toc
+%%
+
+save('symb_abs.mat');
\ No newline at end of file