diff --git a/@SSR/DF_P.m b/@SSR/DF_P.m
index abca5bd..e844e13 100644
--- a/@SSR/DF_P.m
+++ b/@SSR/DF_P.m
@@ -14,5 +14,6 @@
     D = blkdiag(DSC{:});
     BU = blkdiag(BU{:});
 end 
-DF_P = speye( O.nt * O.n ) + BU * get_ConvMtx(O) * D;
+DD = BU * get_ConvMtx(O) * D;
+DF_P = speye( size(DD) ) + DD;
 end
\ No newline at end of file
diff --git a/@SSR/FRC.m b/@SSR/FRC.m
new file mode 100644
index 0000000..0434af1
--- /dev/null
+++ b/@SSR/FRC.m
@@ -0,0 +1,12 @@
+function [om, N] = FRC(O,T_range,ncontsteps)
+% function that computes the FRC of S for a given range
+
+[om,cont_sol,~] = O.sequential_continuation(T_range,'ncontsteps',ncontsteps);
+
+N = nan(length(om),1);
+for j = 1:length(om)
+    if ~isempty(cont_sol{j})
+        N(j) = O.dt*norm(cont_sol{j},'fro');
+    end
+end    
+end
\ No newline at end of file
diff --git a/@SSR/SSM2.m b/@SSR/SSM2.m
new file mode 100644
index 0000000..be074bd
--- /dev/null
+++ b/@SSR/SSM2.m
@@ -0,0 +1,53 @@
+function W = SSM2(O,S)
+% A faster computation of the leading order term of the SSM
+
+n = length(O.M);
+m = O.n;
+W = cell(1,n-m);
+
+dim = (2*m)^2;
+
+A = [zeros(m),eye(m);-O.K1,-O.C1];  % linear coefficient of master system
+B0 = zeros(dim);
+
+delta = eye(2*m);                   % kronecker delta 
+
+% assembly of the 2D equivalent to tensor B as a part of vectorizing the
+% problem (first half that is independent of enslaved mode)
+for i = 1:dim
+    for j = 1:dim
+        q = floor((j-1)/2/m)+1;
+        t = floor((i-1)/2/m)+1;
+        r = j - (q-1) * 2 * m;
+        s = i - (t-1) * 2 * m;
+        B0(i,j) = 2 * A(r,s) * A(q,t) +...
+            A(q,:) * A(:,t) * delta(r,s) +...
+            A(r,:) * A(:,s) * delta(q,t);   % eq. (2.14) in thesis 
+    end
+end
+
+
+% Computing SSM coefficient Wk for each enslaved mode
+for k = 1:n-m
+    B1 = zeros(dim);
+    Sk = zeros(2*m);
+    Sk(1:m,1:m) = S{k};     % this is R_k, nonlinear quadratic coeffs in original eq.
+    
+    % Computing second half of B
+    for i = 1:dim
+        for j = 1:dim
+            q = floor((j-1)/2/m)+1;
+            t = floor((i-1)/2/m)+1;
+            r = j - (q-1) * 2 * m;
+            s = i - (t-1) * 2 * m;
+            B1(i,j) = 2 * (O.zeta2(k)) * (O.omega2(k)) *...
+                (A(q,t) * delta(r,s) + A(r,s) * delta(q,t));
+        end
+    end    
+    
+    B = B0 + B1 + (O.omega2(k))^2 * eye(dim);
+    Wk = -B\Sk(:);                  % quadratic SSM coeffs in vector form
+    W{k} = reshape(Wk,size(Sk));    % converting SSM coeffs back to matrix form
+end
+
+end
\ No newline at end of file
diff --git a/@SSR/SSMcomps.m b/@SSR/SSMcomps.m
new file mode 100644
index 0000000..ea4b408
--- /dev/null
+++ b/@SSR/SSMcomps.m
@@ -0,0 +1,15 @@
+function SSMcomps(O)
+% computes additional material required for SSM computation
+[V1, dd] = eig(full(O.K),full(O.M),'vector');
+[~, ind] = sort(dd);
+V1 = V1(:,ind);
+V2 = V1(:,setdiff(1:length(O.M),O.mode_choice));
+mu2 = diag(V2.' * O.M * V2);
+U2 = V2 * diag( 1./ sqrt(mu2) );
+omega2 = sqrt(diag((U2.' * O.K * U2)));
+zeta2 = diag((U2.' * O.C * U2))./ (2*omega2);
+O.U2store = U2;
+O.omega2store = omega2;
+O.zeta2store = zeta2;
+O.SSMdone = 1;
+end
\ No newline at end of file
diff --git a/@SSR/SSR.m b/@SSR/SSR.m
index 362b5e4..ac9c84f 100644
--- a/@SSR/SSR.m
+++ b/@SSR/SSR.m
@@ -1,6 +1,15 @@
 classdef SSR < handle
-    %SSR Summary of this class goes here
-    %   Detailed explanation goes here
+    %SSR (SteadyStateResponse) class implements the integral equations
+    %approach to the computation of the steady-state response in
+    %(quasi)periodically-forced dynamical systems.
+    
+    % The theory of computation (without model reduction) is decribed in
+    % the following article:
+    %
+    
+    % The theory for model reduction is described in the following article:
+    % G. Buza, S. Jain, G. Haller, Using Spectral Submanifolds for Optimal
+    % Mode Selection in Model Reduction, (2020) Preprint available on arXiv.org
     
     properties
         M               % Mass matrix
@@ -35,6 +44,8 @@
         DLvec           % Derivative dLdT
         Jvec            % Green's function for velocities (calculated over -T:dt:T)
         ConvMtx         % convolution matrix for displacements
+        A               % reformulated convolution matrix
+        Lint            % precomputed integrals of Green's functions
         isupdated       % check if different quantitities are updated according to T
         weights         % Weights vector for higher order integration
         nt              % number of nodes in the time mesh
@@ -51,30 +62,80 @@
         Evinv           % The inverse Fourier transform exponential matrix
         Ekron           % The Fourier transform exponential matrix kronecker product with Identity
         Einvkron        % The inverseFourier transform exponential matrix kronecker product with Identity
+        mode_choice     % Selection of modes specified by the user
+        
+        % required only for the SSM computation in the mode selection
+        % procedure
+        zeta            % The damping coefficients of the master modes
+        omega2          % The eigenfrequencies of the enslaved modes
+        zeta2           % The damping coefficients of the enslaved modes
+        K1              % Stiffness matrix of the enslaved modes (in modal coordinates)
+        C1              % Damping matrix of the enslaved modes (in modal coordinates)
+        U2              % Enslaved eigenmodes
+        % storage for some of these
+        U2store
+        zeta2store
+        omega2store
+        SSMdone         % Boolean that ensures that the eigenvalues are only computed once
     end
     
     methods
-        function O = SSR(M,C,K)
+        function O = SSR(M,C,K,varargin)
             % Basic system properties
             O.M = M;
             O.C = C;
-            O.K = K;
-            O.n = length(M);
+            O.K = K;            
+            if nargin == 3
+                n = length(M);
+                O.mode_choice = 1:n;
+                [VV, dd] = eig(full(K),full(M));
+            else
+                O.mode_choice = varargin{1};
+                [VV, dd] = eigs(K,M,max(O.mode_choice),'SM');
+            end            
+            O.n = length(O.mode_choice);           
             
-            [V, ~] = eig(full(K),full(M));
+            dd = diag(dd);
+            [~, ind] = sort(dd);
+            VV = VV(:,ind);
+            V = VV(:,O.mode_choice);
             mu = diag(V.' * M * V);
             O.U = V * diag( 1./ sqrt(mu) ); % mass normalized VMs
             O.omega0 = sqrt(diag((O.U.' * K * O.U)));
-            zeta = diag((O.U.' * C * O.U))./ (2*O.omega0);
-            O.o = find(zeta>1);
-            O.c = find(zeta==1);
-            O.u = find(zeta<1);
-            lambda1 = (-zeta + sqrt(zeta.^2 - 1)).* O.omega0;
-            lambda2 = (-zeta - sqrt(zeta.^2 - 1)).* O.omega0;
+            O.K1 = O.U.' * K * O.U;
+            O.zeta = diag((O.U.' * C * O.U))./ (2*O.omega0);
+            O.C1 = O.U.' * C * O.U;
+            O.o = find(O.zeta>1);
+            O.c = find(O.zeta==1);
+            O.u = find(O.zeta<1);
+            lambda1 = (-O.zeta + sqrt(O.zeta.^2 - 1)).* O.omega0;
+            lambda2 = (-O.zeta - sqrt(O.zeta.^2 - 1)).* O.omega0;
             O.alpha = real(lambda2);
             O.omega = abs(imag(lambda2));
             O.beta = lambda1;
             O.gamma = lambda2;
+            O.SSMdone = 0;
+        end  
+        
+        function U2 = get.U2(O)
+            if ~O.SSMdone
+                SSMcomps(O);
+            end
+            U2 = O.U2store;
+        end
+        
+        function omega2 = get.omega2(O)
+            if ~O.SSMdone
+                SSMcomps(O);
+            end
+            omega2 = O.omega2store;
+        end
+        
+        function zeta2 = get.zeta2(O)
+            if ~O.SSMdone
+                SSMcomps(O);   
+            end
+            zeta2 = O.zeta2store;
         end
         
         function set.T(O,T)
@@ -188,7 +249,7 @@
                 Jvec = O.Jvec;
             end
         end
-        
+               
         update_Jvec(O)
         
         J = J(O,t,T)
@@ -224,7 +285,7 @@
         
         DF_Q = DF_Q(O,x)
         
-        [x,xd] = LinearResponse(O)
+        [x, xd] = LinearResponse(O)
         
         [x, xd] = Picard(O,varargin)
         
@@ -233,6 +294,12 @@
         [x0,tol,maxiter] = parse_iteration_inputs(O,inputs)
         
         [OMEGA, SOL, PICARD] = sequential_continuation(O,range,varargin)
+               
+        W = SSM2(O,S)
+        
+        [U2,omega2,zeta2] = SSMcomps(O)
+        
+        [om, N] = FRC(O,T_range,ncontsteps)
     end
     
     methods(Static)
diff --git a/examples/ModifiedShawPierrePeriodic.m b/examples/ModifiedShawPierrePeriodic.m
index d3c38c1..b4fa3f7 100644
--- a/examples/ModifiedShawPierrePeriodic.m
+++ b/examples/ModifiedShawPierrePeriodic.m
@@ -2,6 +2,13 @@
 % m*x1dd + k(x1-x2) + c(x1d-x2d) + g*x1^3 = f1(t)
 % m*x2dd + k(2*x2-x1) + c(2*x2d - x1d) = f2(t)
 
+% This example is described in Section 4.1 (Figures 1,2) of the following 
+% article 
+% S. Jain, T. Breunung, G. Haller, Fast Computation of Steady-State 
+% Response for Nonlinear Vibrations of High-Degree-of-Freedom Systems, 
+% Nonlinear Dyn (2019) 97: 313. https://doi.org/10.1007/s11071-019-04971-1
+
+
 clear all
 clc; close all
 
diff --git a/examples/ModifiedShawPierreQuasiPeriodic.m b/examples/ModifiedShawPierreQuasiPeriodic.m
index a2346dd..02703d5 100644
--- a/examples/ModifiedShawPierreQuasiPeriodic.m
+++ b/examples/ModifiedShawPierreQuasiPeriodic.m
@@ -2,6 +2,13 @@
 % m*x1dd + k(x1-x2) + c(x1d-x2d) + g*x1^3 = f1(t)
 % m*x2dd + k(2*x2-x1) + c(2*x2d - x1d) = f2(t)
 
+% This example is described in Section 4.1.2 (Figure 4) of the following 
+% article 
+% S. Jain, T. Breunung, G. Haller, Fast Computation of Steady-State 
+% Response for Nonlinear Vibrations of High-Degree-of-Freedom Systems, 
+% Nonlinear Dyn (2019) 97: 313. https://doi.org/10.1007/s11071-019-04971-1
+
+
 clear all
 clc; close all
 
diff --git a/examples/motiv_example.m b/examples/motiv_example.m
new file mode 100644
index 0000000..acad724
--- /dev/null
+++ b/examples/motiv_example.m
@@ -0,0 +1,101 @@
+%% Mass-spring motivational example for mode selection
+
+% This example is described in Section 4.3 (Figure 2) of the following 
+% article 
+% G. Buza, S. Jain, G. Haller, Using Spectral Submanifolds for Optimal 
+% Mode Selection in Model Reduction, (2020) Preprint available on arXiv.org
+
+clear 
+clc
+close all
+
+%% properties
+M = eye(3);
+D1 = 0.01;
+D2 = 0.02;
+D3 = 0.08;
+ome1 = 2;
+ome2 = 3;
+ome3 = 5;
+K = diag([ome1^2 ome2^2 ome3^2]);
+C = 2*diag([ome1*D1 ome2*D2 ome3*D3]);
+
+% Nonlinearity
+S = @(x)SP_nonlinearity3(x,ome1,ome2,ome3);
+DS = @(x)SP_nonlinearity_derv3(x,ome1,ome2,ome3);
+
+% Loading
+A = 0.02;                     % loading amplitude       
+f0 = [1;0;0];                 % loading shape
+f = @(t,T)A*f0*sin(2*pi*t/T); % loading function
+
+%% Mode selection
+
+SS = SSR(M,C,K,1);
+
+% SSM computation
+R = cell(1,2);
+R{1} = ome2^2/2;
+R{2} = ome3^2/2;
+W = SS.SSM2(R);
+disp(['norm(W_2)=' num2str(norm(W{1})) ', norm(W_3)=' num2str(norm(W{2}))])
+[~, ind] = max([norm(W{1}) norm(W{2})]);
+
+% master mode sets considered
+I_1 = [1 2];        % mode selection 1
+I_2 = [1 ind+1];    % mode selection 2
+
+%% SSR Package
+SSfull = SSR(M,C,K,[1 2 3]);   % full system
+SSred1 = SSR(M,C,K,I_1);     % reduced to I_1
+SSred2 = SSR(M,C,K,I_2);     % reduced to I_2
+SSfull.S = S;
+SSfull.DS = DS; 
+SSfull.f = f;
+SSred1 = copyfromfull(SSred1,SSfull);
+SSred2 = copyfromfull(SSred2,SSfull);
+
+%% Sequential continuation
+omega1 = min(SSfull.omega);
+Omega_range = [0.7 1.3]*omega1;
+T_range = 2*pi./Omega_range;
+ncontsteps = 100;
+
+[Omega_arrayf,Nf] = SSfull.FRC(T_range,ncontsteps);
+[Omega_array1,N1] = SSred1.FRC(T_range,ncontsteps);
+[Omega_array2,N2] = SSred2.FRC(T_range,ncontsteps);
+
+Nlin = nan(length(Omega_arrayf),1);
+for j = 1:length(Omega_arrayf)
+    SSfull.T = 2*pi/Omega_arrayf(j);
+    [x_lin, xd_lin] = SSfull.LinearResponse();
+    LinSol = [x_lin; xd_lin];
+    Nlin(j) = SSfull.dt*norm(LinSol,'fro');
+end
+
+figure(2); semilogy(Omega_arrayf,Nf,'-k','DisplayName', 'Full','linewidth',1); axis tight; grid on; hold on;
+xlabel('$$\Omega$$ [rad/s]'); ylabel('$$||q||_2$$')
+legend('show')
+semilogy(Omega_arrayf,Nlin,'-b', 'DisplayName', 'Linear','linewidth',1);
+semilogy(Omega_array1,N1,'--g', 'DisplayName', 'Reduced $$I_1$$','linewidth',2);
+semilogy(Omega_array2,N2,'--r','DisplayName', 'Reduced $$I_2$$','linewidth',2);
+
+
+function f = SP_nonlinearity3(x,o1,o2,o3)
+f = [o1^2/2*(3*x(1)^2+x(2)^2+x(3)^2)+o2^2*x(1)*x(2)+o3^2*x(1)*x(3)+(o1^2+o2^2+o3^2)/2*x(1)*(x(1)^2+x(2)^2+x(3)^2);...
+    o2^2/2*(3*x(2)^2+x(1)^2+x(3)^2)+o1^2*x(1)*x(2)+o3^2*x(2)*x(3)+(o1^2+o2^2+o3^2)/2*x(2)*(x(1)^2+x(2)^2+x(3)^2);...
+    o3^2/2*(3*x(3)^2+x(2)^2+x(1)^2)+o2^2*x(3)*x(2)+o1^2*x(1)*x(3)+(o1^2+o2^2+o3^2)/2*x(3)*(x(1)^2+x(2)^2+x(3)^2) ];
+end
+
+function Df = SP_nonlinearity_derv3(x,o1,o2,o3)
+Df = zeros(3);
+Df(1,1) = 3*o1^2*x(1)+o2^2*x(2)+o3^2*x(3) + (o1^2+o2^2+o3^2)/2*(3*x(1)^2 + x(2)^2+x(3)^2);
+Df(2,2) = 3*o2^2*x(2)+o1^2*x(1)+o3^2*x(3) + (o1^2+o2^2+o3^2)/2*(3*x(2)^2 + x(1)^2+x(3)^2);
+Df(3,3) = (3*o3^2*x(3)+o2^2*x(2)+o1^2*x(1) + (o1^2+o2^2+o3^2)/2*(3*x(3)^2 + x(2)^2+x(1)^2));
+Df(1,2) = o1^2*x(2)+o2^2*x(1)+(o1^2+o2^2+o3^2)*x(1)*x(2);
+Df(1,3) = o1^2*x(3)+o3^2*x(1)+(o1^2+o2^2+o3^2)*x(1)*x(3);
+Df(2,3) = o2^2*x(3)+o3^2*x(2)+(o1^2+o2^2+o3^2)*x(2)*x(3);
+Df(3,2) = Df(2,3);
+Df(2,1) = Df(1,2);
+Df(3,1) = Df(1,3);
+end
\ No newline at end of file
diff --git a/examples/ndof_oscillator_chain/OscillatorChain.m b/examples/ndof_oscillator_chain/OscillatorChain.m
index d1df4fb..d6007fc 100644
--- a/examples/ndof_oscillator_chain/OscillatorChain.m
+++ b/examples/ndof_oscillator_chain/OscillatorChain.m
@@ -1,4 +1,10 @@
 %% Oscillator chain with coupled nonlinearities
+% This example is described in Section 4.2 (Figure 6) of the following 
+% article 
+% S. Jain, T. Breunung, G. Haller, Fast Computation of Steady-State 
+% Response for Nonlinear Vibrations of High-Degree-of-Freedom Systems, 
+% Nonlinear Dyn (2019) 97: 313. https://doi.org/10.1007/s11071-019-04971-1
+
 clear all
 clc
 
diff --git a/examples/nonsmooth/NonSmoothOscillator.m b/examples/nonsmooth/NonSmoothOscillator.m
index c81b7c3..dd6d739 100644
--- a/examples/nonsmooth/NonSmoothOscillator.m
+++ b/examples/nonsmooth/NonSmoothOscillator.m
@@ -1,4 +1,11 @@
 %% 2-DOF oscillator with non smooth nonlinearity
+% This example is described in Section 4.1.3 (Figure 5) of the following 
+% article 
+% S. Jain, T. Breunung, G. Haller, Fast Computation of Steady-State 
+% Response for Nonlinear Vibrations of High-Degree-of-Freedom Systems, 
+% Nonlinear Dyn (2019) 97: 313. https://doi.org/10.1007/s11071-019-04971-1
+
+
 clear all
 clc; close all
 
diff --git a/examples/vonKarman/example_curvedbeam.m b/examples/vonKarman/example_curvedbeam.m
new file mode 100644
index 0000000..ca31e7a
--- /dev/null
+++ b/examples/vonKarman/example_curvedbeam.m
@@ -0,0 +1,99 @@
+%% Mode Selection on Von Karman curved Beam Example
+% This example is described in Section 6.2 (Figure 8) of the following 
+% article. 
+% G. Buza, S. Jain, G. Haller, Using Spectral Submanifolds for Optimal 
+% Mode Selection in Model Reduction, (2020) Preprint available on arXiv.org
+
+% The finite element model is taken from the following article. 
+% Jain, S., Tiso, P., & Haller, G. (2018). Exact nonlinear model reduction 
+% for a von Kármán beam: slow-fast decomposition and spectral submanifolds. 
+% Journal of Sound and Vibration, 423, 195–211. 
+% https://doi.org/10.1016/J.JSV.2018.01.049
+
+clear 
+clc
+close all
+
+%% parameters
+epsilon = 7e-3; % beam thickness-to-length ratio 
+midp_height = 0.005; % height of midpoint relative to the ends (measure of curvature of the beam)
+[Geo] = Geometry(epsilon,midp_height);
+
+nElements = 10;
+BC = 'B'; % boundary condition: simply supported
+
+%% build finite element model
+[Misc,model] = Beam_Model(Geo,nElements,BC); 
+
+n = length(model.freeDOFs); % total number of degrees of freedom
+
+%% Mass Matrix
+M = model.M(model.freeDOFs,model.freeDOFs);
+
+%% Stiffness Matrix
+K = model.K(model.freeDOFs,model.freeDOFs);
+
+%% Damping Matrix
+C = model.C(model.freeDOFs,model.freeDOFs);
+
+%% Nonlinearity
+S =@(x) NonlinearityVK(model,x);
+DS =@(x) NonlinearityJacobianVK(model,x);
+
+A = 80; % amplitude for which the FRC is computed
+f = @(t,T) A*model.loads(model.freeDOFs) * sin(2*pi*t/T); % loading function
+
+%% Steady-State tool
+SSfull = SSR(M,C,K,1:n);        % full system
+
+%% Mode selection
+
+% Input parameters for mode selection
+param.curvtolerance = 0.05;
+param.nmodes = 10;
+param.type = 1; % set to 0 if param.nmodes is the max number of modes,
+% set to 1 if param.nmodes is the desired number of modes
+param.example = 'simple'; % always use simple if enough memory is available to store X below
+X = S11(model); % this conains the quadratic nonlinearities for all dofs (X_i).
+% X_i is symmetric, and can be found from the form M\ddot{q} + C\dot{q} +
+% Kq +\Gamma(q) = F(t) (i denotes the specific row of this equation)
+
+% master mode sets considered
+modesel1 = 1:10; 
+modesel2 = modeselect(SSfull,param,X);
+SSred1 = SSR(M,C,K,modesel1);   % reduced to I_1
+SSred2 = SSR(M,C,K,modesel2);   % reduced to I_2
+
+SSfull.S = S;
+SSfull.DS = DS;
+SSfull.f = f;
+SSred1 = copyfromfull(SSred1,SSfull);
+SSred2 = copyfromfull(SSred2,SSfull);
+
+%% Sequential continuation
+omega1 = SSfull.omega(1);
+Omega_range = [0.7 1.3]*omega1;             
+T_range = 2*pi./Omega_range;
+ncontsteps = 100;
+
+[Omega_arrayf,Nf] = SSfull.FRC(T_range,ncontsteps);
+[Omega_array1,N1] = SSred1.FRC(T_range,ncontsteps);
+[Omega_array2,N2] = SSred2.FRC(T_range,ncontsteps);
+
+Nlin = nan(length(Omega_arrayf),1);
+for j = 1:length(Omega_arrayf)
+    SSfull.T = 2*pi/Omega_arrayf(j);
+    [x_lin, xd_lin] = SSfull.LinearResponse();
+    LinSol = [x_lin; xd_lin];
+    Nlin(j) = SSfull.dt*norm(LinSol,'fro');
+end
+
+figure(2); semilogy(Omega_arrayf,Nf,'-k','DisplayName', 'Full','linewidth',1); axis tight; grid on; hold on;
+xlabel('$$\Omega$$ [rad/s]'); ylabel('$$||q||_2$$')
+legend('show')
+semilogy(Omega_arrayf,Nlin,'-b', 'DisplayName', 'Linear','linewidth',1);
+semilogy(Omega_array1,N1,'--g', 'DisplayName', 'Reduced $$I_1$$','linewidth',2);
+semilogy(Omega_array2,N2,'--r','DisplayName', 'Reduced $$I_2$$','linewidth',2);
+
+
+
diff --git a/examples/vonKarman/example_straightbeam.m b/examples/vonKarman/example_straightbeam.m
new file mode 100644
index 0000000..25c298f
--- /dev/null
+++ b/examples/vonKarman/example_straightbeam.m
@@ -0,0 +1,136 @@
+%% Mode selection on Von Karman Straight Beam example
+% This example is described in Section 6.1 (Figure 5) of the following 
+% article. 
+% G. Buza, S. Jain, G. Haller, Using Spectral Submanifolds for Optimal 
+% Mode Selection in Model Reduction, (2020) Preprint available on arXiv.org
+
+% The finite element model is taken from the following article. 
+% Jain, S., Tiso, P., & Haller, G. (2018). Exact nonlinear model reduction 
+% for a von Kármán beam: slow-fast decomposition and spectral submanifolds. 
+% Journal of Sound and Vibration, 423, 195–211. 
+% https://doi.org/10.1016/J.JSV.2018.01.049
+
+clear 
+clc
+close all
+%% parameters
+epsilon = 1e-3; % beam thickness-to-length ratio 
+midp_height = 0; % height of midpoint relative to the ends (measure of curvature of the beam)
+[Geo] = Geometry(epsilon,midp_height);
+
+nElements = 10;
+BC = 'B'; % boundary condition: simply supported
+
+%% build finite element model
+[Misc,model] = Beam_Model(Geo,nElements,BC); 
+
+n = length(model.freeDOFs); % total number of degrees of freedom
+plot_dof = 9; % the degree of freedom at which the response is plotted
+%% Mass Matrix
+M = model.M(model.freeDOFs,model.freeDOFs);
+
+%% Stiffness Matrix
+K = model.K(model.freeDOFs,model.freeDOFs);
+
+%% Damping Matrix
+C = model.C(model.freeDOFs,model.freeDOFs);
+
+%% Steady-State tool
+SS = SSR(M,C,K,1:n);   % full system
+
+%% Mode selection
+
+% Input parameters for mode selection
+param.curvtolerance = 0.05;
+param.nmodes = 10;
+param.type = 0; % set to 0 if param.nmodes is the max number of modes,
+                % set to 1 if param.nmodes is the desired number of modes
+param.example = 'simple'; % always use simple if enough memory is available to store X below
+param.initialmodes = [1 2 3 4 5]; % optional input modes
+X = S11(model); % this conains the quadratic nonlinearities for all dofs (X_i).
+                % X_i is symmetric, and can be found from the form M\ddot{q} + C\dot{q} +
+                % Kq +\Gamma(q) = F(t) (i denotes the specific row of this equation)
+
+% master mode sets considered
+I_1 = 1:10; 
+I_2 = modeselect(SS,param,X);
+SS1 = SSR(M,C,K,I_1);  % reduced to I_1   
+SS2 = SSR(M,C,K,I_2);  % reduced to I_2
+
+%% Incrementing amplitude
+T = 2*pi/26;
+A = [10:5:35 40:4:172 174 176:1:200 200:1:230]* 10^(-2); % amplitude vector
+SS.S = @(x) NonlinearityVK(model,x);
+SS.DS = @(x) NonlinearityJacobianVK(model,x); 
+
+tol = 1e-6; % Newton Raphson error tolerance
+maxiter = 60;
+N2 = nan(length(A),1);
+N1 = nan(length(A),1);
+NF = nan(length(A),1);
+NL = nan(length(A),1);
+
+for j = 1:length(A)
+    Fext = @(t,T) A(j)*model.loads(model.freeDOFs) * sin(2*pi*t/T);
+    SS.f = Fext;          
+    SS.T = T;  
+    SS.n_steps = 100;
+    SS1 = copyfromfull(SS1,SS);
+    SS2 = copyfromfull(SS2,SS);
+
+    if j == 1
+        [x2,xd2] = SS1.NewtonRaphson('tol',tol,'maxiter',maxiter); % find out response using Picard iteration
+    else
+        [x2,xd2] = SS1.NewtonRaphson('init',x2,'tol',tol,'maxiter',maxiter);
+    end
+    N1(j) = SS1.dt*norm([x2; xd2],'fro');
+
+    if j == 1
+        [x3,xd3] = SS2.NewtonRaphson('tol',tol,'maxiter',maxiter); % find out response using Picard iteration
+    else
+        [x3,xd3] = SS2.NewtonRaphson('init',x3,'tol',tol,'maxiter',maxiter);
+    end
+    N2(j) = SS2.dt*norm([x3; xd3],'fro');
+
+
+    if j == 1
+        [xf,xdf] = SS.NewtonRaphson('tol',tol,'maxiter',maxiter);
+    else
+        [xf,xdf] = SS.NewtonRaphson('init',xf,'tol',tol,'maxiter',maxiter);
+    end
+
+    NF(j) = SS.dt*norm([xf; xdf],'fro');
+
+    [x_lin, xd_lin] = SS.LinearResponse();
+    NL(j) = SS.dt*norm([x_lin; xd_lin],'fro');
+    
+    disp(vpa(A(j)))
+end
+
+figure(3);
+
+%% final amplitude time history
+plot(SS.t, xf(plot_dof,:), 'k', 'DisplayName', 'Full','linewidth',1);
+axis tight;
+legend('show')
+ylim([-3.1*10^(-3) 3.1*10^(-3)])
+xlabel('time (seconds)');
+ylabel('displacement $$w_{4}$$ [m]')
+hold on;
+plot(SS.t, x_lin(plot_dof,:), 'b', 'DisplayName', 'Linear','linewidth',1);
+plot(SS1.t, x2(plot_dof,:),'--g', 'DisplayName', 'Reduced $$I_1$$','linewidth',2);
+plot(SS2.t, x3(plot_dof,:),'--r', 'DisplayName', 'Reduced $$I_2$$','linewidth',2);
+
+%% forcing amp - response amp curve
+figure
+plot(A,NF,'-k','DisplayName', 'Full','linewidth',1); axis tight; hold on;
+xlabel('$$F$$ [N]'); ylabel('$$||q||_2$$')
+lgd = legend('show');
+lgd.Location = 'southeast';
+plot(A,NL,'-b', 'DisplayName', 'Linear','linewidth',1);
+plot(A,N1,'--g', 'DisplayName', 'Reduced $$I_1$$','linewidth',2);
+plot(A,N2,'--r','DisplayName', 'Reduced $$I_2$$','linewidth',2);
+ylim([0 0.032])
+
+
+
diff --git a/examples/vonKarman/model/AssembleLinearStiffness.m b/examples/vonKarman/model/AssembleLinearStiffness.m
new file mode 100644
index 0000000..d27a5b0
--- /dev/null
+++ b/examples/vonKarman/model/AssembleLinearStiffness.m
@@ -0,0 +1,33 @@
+function [K] = AssembleLinearStiffness(model)
+% The finite element model is taken from the following article: 
+% Jain, S., Tiso, P., & Haller, G. (2018). Exact nonlinear model reduction 
+% for a von Kármán beam: slow-fast decomposition and spectral submanifolds. 
+% Journal of Sound and Vibration, 423, 195–211. 
+% https://doi.org/10.1016/J.JSV.2018.01.049
+
+ndof = model.nDOF;
+nelem = model.nElements;
+I = zeros(nelem*6*6,1);
+J = zeros(nelem*6*6,1);
+KK = zeros(nelem*6*6,1);
+
+el = 0;
+% F = zeros(ndof,1);
+for i = 1:model.nElements
+    index = 3*(i-1)+1:3*(i+1);
+    
+    Te = compute_rotation_matrix(model.nodes(i:i+1,:));
+
+    [Kel,~] = ElementStiffness(model,zeros(6,1));
+    Kel = Te.' * Kel * Te; 
+           
+%      F(index) = F(index)+ Fel;
+     for col = 1:6
+        I(el+6*(col-1)+1:el+6*col,1) = ones(6,1)*index(col);
+        J(el+6*(col-1)+1:el+6*col,1) = index;
+        KK(el+6*(col-1)+1:el+6*col,1) = Kel(col,:);
+    end
+    el = el + 36; 
+  
+end
+K = sparse(I,J,KK,ndof,ndof);
diff --git a/examples/vonKarman/model/AssembleMass.m b/examples/vonKarman/model/AssembleMass.m
new file mode 100644
index 0000000..607fbd7
--- /dev/null
+++ b/examples/vonKarman/model/AssembleMass.m
@@ -0,0 +1,30 @@
+function M = AssembleMass(model)
+% The finite element model is taken from the following article: 
+% Jain, S., Tiso, P., & Haller, G. (2018). Exact nonlinear model reduction 
+% for a von Kármán beam: slow-fast decomposition and spectral submanifolds. 
+% Journal of Sound and Vibration, 423, 195–211. 
+% https://doi.org/10.1016/J.JSV.2018.01.049
+
+ndof = model.nDOF;
+nelem = model.nElements;
+I = zeros(nelem*6*6,1);
+J = zeros(nelem*6*6,1);
+KK = zeros(nelem*6*6,1);
+el = 0;
+for i = 1:model.nElements
+    index = 3*(i-1)+1:3*(i+1);
+    
+    Te = compute_rotation_matrix(model.nodes(i:i+1,:));
+    
+    Mel = ElementMass(model);
+    Mel = Te.' * Mel * Te; 
+        
+    for col = 1:6
+        I(el+6*(col-1)+1:el+6*col,1) = ones(6,1)*index(col);
+        J(el+6*(col-1)+1:el+6*col,1) = index;
+        KK(el+6*(col-1)+1:el+6*col,1) = Mel(col,:);
+    end
+    el = el + 36; 
+  
+end
+M = sparse(I,J,KK,ndof,ndof);
diff --git a/examples/vonKarman/model/AssembleTangentStiffness.m b/examples/vonKarman/model/AssembleTangentStiffness.m
new file mode 100644
index 0000000..dcd42e3
--- /dev/null
+++ b/examples/vonKarman/model/AssembleTangentStiffness.m
@@ -0,0 +1,35 @@
+function [K,F] = AssembleTangentStiffness(model,u)
+% The finite element model is taken from the following article: 
+% Jain, S., Tiso, P., & Haller, G. (2018). Exact nonlinear model reduction 
+% for a von Kármán beam: slow-fast decomposition and spectral submanifolds. 
+% Journal of Sound and Vibration, 423, 195–211. 
+% https://doi.org/10.1016/J.JSV.2018.01.049
+
+ndof = model.nDOF;
+nelem = model.nElements;
+I = zeros(nelem*6*6,1);
+J = zeros(nelem*6*6,1);
+KK = zeros(nelem*6*6,1);
+
+el = 0;
+F = zeros(ndof,1);
+for i = 1:model.nElements
+    index = 3*(i-1)+1:3*(i+1);
+    
+    Te = compute_rotation_matrix(model.nodes(i:i+1,:));
+    
+    qe = Te * u(index); 
+    [Kel,Fel] = ElementStiffness(model,qe);
+    
+    Kel = Te.' * Kel * Te; 
+    F(index) = F(index)+ Te.' * Fel;
+    
+    for col = 1:6
+        I(el+6*(col-1)+1:el+6*col,1) = ones(6,1)*index(col);
+        J(el+6*(col-1)+1:el+6*col,1) = index;
+        KK(el+6*(col-1)+1:el+6*col,1) = Kel(col,:);
+    end
+    el = el + 36; 
+  
+end
+K = sparse(I,J,KK,ndof,ndof);
diff --git a/examples/vonKarman/model/Beam_Model.m b/examples/vonKarman/model/Beam_Model.m
new file mode 100644
index 0000000..3ca72f9
--- /dev/null
+++ b/examples/vonKarman/model/Beam_Model.m
@@ -0,0 +1,89 @@
+function [Misc,model] = Beam_Model(Geo,nElements,B)
+% This finite element model is taken from the following article: 
+% Jain, S., Tiso, P., & Haller, G. (2018). Exact nonlinear model reduction 
+% for a von Kármán beam: slow-fast decomposition and spectral submanifolds. 
+% Journal of Sound and Vibration, 423, 195–211. 
+% https://doi.org/10.1016/J.JSV.2018.01.049
+
+Damping = 2;    % Structural damping : 0 - Nil,  1 - Rayleigh, 2 - Linear viscoelastic
+
+% Assumed a curved beam with N equidistant elements along the x axis
+model.nElements = nElements;
+model.nNodes = nElements +1;
+model.nDOF = model.nNodes*3;
+model.dx = Geo.l/model.nElements;
+model.b = Geo.b;
+
+model.nodes = [(0:model.dx:Geo.l).' zeros(model.nNodes,1)];
+if Geo.a~=0
+    R = 1/2*(Geo.a^2+(Geo.l/2)^2)/Geo.a;
+    theta0 = asin(Geo.l/2/R);
+    
+    for i = 1:size(model.nodes,1)
+        th = asin((model.nodes(i,1)-Geo.l/2)/R);
+        model.nodes(i,2) = R *  cos(th) - R * cos(theta0);
+    end
+end
+
+% Material Properties (Aluminium)
+model.E = 70e9;     % Youngs Modulus (Pa)
+model.A = Geo.A;    % Cross section area of Beam
+model.I = Geo.I;    % Moment of inertia of beam
+model.rho = 2700;   % Density
+model.kappa = 1e8;  % Material damping for Kelvin Voigt model
+     
+
+% Loading
+model.pw = 1;      % loading factor along the transverse direction
+model.pu = 0;      % loading factor along the axial direction
+
+model.loads = sparse(model.nDOF,1);
+l = ones(model.nNodes,1);
+l(2:end-1) = 2;
+model.uDOF = 1:3:model.nDOF; % axial degrees of freedom
+model.wDOF = setdiff(1:model.nDOF,model.uDOF); % transverse DOFs 
+model.loads(2:3:end) = model.pw*(model.dx/2)*l; % 
+model.loads(model.uDOF) = model.pu*(model.dx/2)*l;
+
+% Boundary conditions
+model.fixDOFs = BoundaryCoundition(model,B);
+model.freeDOFs = setdiff(1:model.nDOF, model.fixDOFs);
+model.uDOF = setdiff(model.uDOF,model.fixDOFs);
+model.wDOF = setdiff(model.wDOF,model.fixDOFs);
+
+% Stiffness
+model.K = AssembleLinearStiffness(model);
+model.M = AssembleMass(model);
+
+%
+model.modes = 1:5;
+model.m = length(model.modes);
+% Vibration problem
+Mred = model.M(model.freeDOFs,model.freeDOFs);
+model.Ks = model.K;
+Kred = model.K(model.freeDOFs,model.freeDOFs);
+[model.V,D] = eigs(Kred,Mred,max(model.modes),'SM');   % Return M smallest eigenvalues and vectors
+[omega2,V] = sortVMsVK(model.V,Mred,D);
+disp('First 5 undamped Natural Frequencies [Hz]')
+disp(sqrt(omega2)/(2*pi))
+model.omega2 = omega2(model.modes);
+model.VMs = zeros(model.nDOF,model.m);
+model.VMs(model.freeDOFs,:) = V(:,model.modes);
+Misc.omega = sqrt(model.omega2);
+
+
+figure
+PlotDeformation(model,model.VMs(:,1),0.1);
+
+switch Damping
+    case 1
+        % Rayleigh Damping
+        W =   sqrt(model.omega2(1:2));
+        a = [W(1) 1/W(1);W(2) 1/W(2)]\[0.02;0.02];
+        model.C = (a(2)*model.M + a(1)*model.K);
+               
+    case 2
+        % ViscoElastic Damping with linear strain rate
+        model.C = (model.kappa/model.E)*model.Ks;        
+end
+
diff --git a/examples/vonKarman/model/BoundaryCoundition.m b/examples/vonKarman/model/BoundaryCoundition.m
new file mode 100644
index 0000000..4fa2e49
--- /dev/null
+++ b/examples/vonKarman/model/BoundaryCoundition.m
@@ -0,0 +1,11 @@
+function [fixDOFs] = BoundaryCoundition(model,S)
+
+switch S
+    case 'C'    % Cantilever
+        fixDOFs = 1:3;
+    case 'B'    % Simply supported
+        fixDOFs = [1 2 model.nDOF-2 model.nDOF-1];
+    case 'DC'  % Doubly clamped
+        fixDOFs = [1:3 (model.nDOF-2):model.nDOF];
+end
+        
\ No newline at end of file
diff --git a/examples/vonKarman/model/ElementMass.m b/examples/vonKarman/model/ElementMass.m
new file mode 100644
index 0000000..ce6a432
--- /dev/null
+++ b/examples/vonKarman/model/ElementMass.m
@@ -0,0 +1,22 @@
+function M = ElementMass(prop)
+% The finite element model is taken from the following article: 
+% Jain, S., Tiso, P., & Haller, G. (2018). Exact nonlinear model reduction 
+% for a von Kármán beam: slow-fast decomposition and spectral submanifolds. 
+% Journal of Sound and Vibration, 423, 195–211. 
+% https://doi.org/10.1016/J.JSV.2018.01.049
+
+rho = prop.rho;
+A   = prop.A;
+l   = prop.dx;
+
+
+M = sparse(6,6);
+
+M([1 4],[1 4]) = ((rho*A*l)/6)*[2 1;
+                                        1 2];
+                              
+M([2 3 5 6],[2 3 5 6]) = ((rho*A*l)/420)*[156      22*l       54      -13*l;
+                                      22*l     4*l^2      13*l    -3*l^2;
+                                      54       13*l       156     -22*l;  
+                                      -13*l     -3*l^2     -22*l   4*l^2];
+
diff --git a/examples/vonKarman/model/ElementStiffness.m b/examples/vonKarman/model/ElementStiffness.m
new file mode 100644
index 0000000..b6ede00
--- /dev/null
+++ b/examples/vonKarman/model/ElementStiffness.m
@@ -0,0 +1,34 @@
+function [K,F] = ElementStiffness(model,q)
+% The finite element model is taken from the following article: 
+% Jain, S., Tiso, P., & Haller, G. (2018). Exact nonlinear model reduction 
+% for a von Kármán beam: slow-fast decomposition and spectral submanifolds. 
+% Journal of Sound and Vibration, 423, 195–211. 
+% https://doi.org/10.1016/J.JSV.2018.01.049
+
+E = model.E;
+A = model.A;
+dx = model.dx;
+I = model.I;
+
+u1 = q(1);
+w1 = q(2);
+t1 = q(3);
+u2 = q(4);
+w2 = q(5);
+t2 = q(6);
+
+
+K = [                                                                 (A*E)/dx,                                                                                                                                                                                                                                        -((A*E*(36*w1 - 36*w2))/30 + (A*E*dx*(3*t1 + 3*t2))/30)/dx^2,                                                                                                                                                                                                                   - (A*E*(4*t1 - t2))/30 - (A*E*(3*w1 - 3*w2))/(30*dx),                                                                        -(A*E)/dx,                                                                                                                                                                                                                                         ((A*E*(36*w1 - 36*w2))/30 + (A*E*dx*(3*t1 + 3*t2))/30)/dx^2,                                                                                                                                                                                                                     (A*E*(t1 - 4*t2))/30 - (A*E*(3*w1 - 3*w2))/(30*dx);
+    (6*A*E*w2)/(5*dx^2) - (A*E*t2)/(10*dx) - (6*A*E*w1)/(5*dx^2) - (A*E*t1)/(10*dx), (12*E*I)/dx^3 - (6*A*E*u1)/(5*dx^2) + (6*A*E*u2)/(5*dx^2) + (9*A*E*t1^2)/(70*dx) + (9*A*E*t2^2)/(70*dx) + (108*A*E*w1^2)/(35*dx^3) + (108*A*E*w2^2)/(35*dx^3) + (27*A*E*t1*w1)/(35*dx^2) - (27*A*E*t1*w2)/(35*dx^2) + (27*A*E*t2*w1)/(35*dx^2) - (27*A*E*t2*w2)/(35*dx^2) - (216*A*E*w1*w2)/(35*dx^3),                           (6*E*I)/dx^2 - (3*A*E*t1^2)/280 + (3*A*E*t2^2)/280 - (A*E*u1)/(10*dx) + (A*E*u2)/(10*dx) + (27*A*E*w1^2)/(70*dx^2) + (27*A*E*w2^2)/(70*dx^2) + (3*A*E*t1*t2)/140 + (9*A*E*t1*w1)/(35*dx) - (9*A*E*t1*w2)/(35*dx) - (27*A*E*w1*w2)/(35*dx^2), (A*E*t1)/(10*dx) + (A*E*t2)/(10*dx) + (6*A*E*w1)/(5*dx^2) - (6*A*E*w2)/(5*dx^2), (6*A*E*u1)/(5*dx^2) - (12*E*I)/dx^3 - (6*A*E*u2)/(5*dx^2) - (9*A*E*t1^2)/(70*dx) - (9*A*E*t2^2)/(70*dx) - (108*A*E*w1^2)/(35*dx^3) - (108*A*E*w2^2)/(35*dx^3) - (27*A*E*t1*w1)/(35*dx^2) + (27*A*E*t1*w2)/(35*dx^2) - (27*A*E*t2*w1)/(35*dx^2) + (27*A*E*t2*w2)/(35*dx^2) + (216*A*E*w1*w2)/(35*dx^3),                           (6*E*I)/dx^2 + (3*A*E*t1^2)/280 - (3*A*E*t2^2)/280 - (A*E*u1)/(10*dx) + (A*E*u2)/(10*dx) + (27*A*E*w1^2)/(70*dx^2) + (27*A*E*w2^2)/(70*dx^2) + (3*A*E*t1*t2)/140 + (9*A*E*t2*w1)/(35*dx) - (9*A*E*t2*w2)/(35*dx) - (27*A*E*w1*w2)/(35*dx^2);
+    (A*E*t2)/30 - (2*A*E*t1)/15 - (A*E*w1)/(10*dx) + (A*E*w2)/(10*dx),                                                       (6*E*I)/dx^2 - (3*A*E*t1^2)/280 + (3*A*E*t2^2)/280 - (A*E*u1)/(10*dx) + (A*E*u2)/(10*dx) + (27*A*E*w1^2)/(70*dx^2) + (27*A*E*w2^2)/(70*dx^2) + (3*A*E*t1*t2)/140 + (9*A*E*t1*w1)/(35*dx) - (9*A*E*t1*w2)/(35*dx) - (27*A*E*w1*w2)/(35*dx^2), (3*A*E*dx*t1^2)/35 - (3*A*E*dx*t1*t2)/140 - (3*A*E*t1*w1)/140 + (3*A*E*t1*w2)/140 + (A*E*dx*t2^2)/140 + (3*A*E*t2*w1)/140 - (3*A*E*t2*w2)/140 + (9*A*E*w1^2)/(70*dx) - (9*A*E*w1*w2)/(35*dx) + (9*A*E*w2^2)/(70*dx) - (2*A*E*u1)/15 + (2*A*E*u2)/15 + (4*E*I)/dx,             (2*A*E*t1)/15 - (A*E*t2)/30 + (A*E*w1)/(10*dx) - (A*E*w2)/(10*dx),                                                       (3*A*E*t1^2)/280 - (6*E*I)/dx^2 - (3*A*E*t2^2)/280 + (A*E*u1)/(10*dx) - (A*E*u2)/(10*dx) - (27*A*E*w1^2)/(70*dx^2) - (27*A*E*w2^2)/(70*dx^2) - (3*A*E*t1*t2)/140 - (9*A*E*t1*w1)/(35*dx) + (9*A*E*t1*w2)/(35*dx) + (27*A*E*w1*w2)/(35*dx^2),                                                                        (A*E*u1)/30 - (A*E*u2)/30 + (2*E*I)/dx - (3*A*E*dx*t1^2)/280 - (3*A*E*dx*t2^2)/280 + (3*A*E*t1*w1)/140 - (3*A*E*t1*w2)/140 + (3*A*E*t2*w1)/140 - (3*A*E*t2*w2)/140 + (A*E*dx*t1*t2)/70;
+    -(A*E)/dx,                                                                                                                                                                                                                                         ((A*E*(36*w1 - 36*w2))/30 + (A*E*dx*(3*t1 + 3*t2))/30)/dx^2,                                                                                                                                                                                                                     (A*E*(4*t1 - t2))/30 + (A*E*(3*w1 - 3*w2))/(30*dx),                                                                         (A*E)/dx,                                                                                                                                                                                                                                        -((A*E*(36*w1 - 36*w2))/30 + (A*E*dx*(3*t1 + 3*t2))/30)/dx^2,                                                                                                                                                                                                                     (A*E*(3*w1 - 3*w2))/(30*dx) - (A*E*(t1 - 4*t2))/30;
+    (A*E*t1)/(10*dx) + (A*E*t2)/(10*dx) + (6*A*E*w1)/(5*dx^2) - (6*A*E*w2)/(5*dx^2), (6*A*E*u1)/(5*dx^2) - (12*E*I)/dx^3 - (6*A*E*u2)/(5*dx^2) - (9*A*E*t1^2)/(70*dx) - (9*A*E*t2^2)/(70*dx) - (108*A*E*w1^2)/(35*dx^3) - (108*A*E*w2^2)/(35*dx^3) - (27*A*E*t1*w1)/(35*dx^2) + (27*A*E*t1*w2)/(35*dx^2) - (27*A*E*t2*w1)/(35*dx^2) + (27*A*E*t2*w2)/(35*dx^2) + (216*A*E*w1*w2)/(35*dx^3),                           (3*A*E*t1^2)/280 - (6*E*I)/dx^2 - (3*A*E*t2^2)/280 + (A*E*u1)/(10*dx) - (A*E*u2)/(10*dx) - (27*A*E*w1^2)/(70*dx^2) - (27*A*E*w2^2)/(70*dx^2) - (3*A*E*t1*t2)/140 - (9*A*E*t1*w1)/(35*dx) + (9*A*E*t1*w2)/(35*dx) + (27*A*E*w1*w2)/(35*dx^2), (6*A*E*w2)/(5*dx^2) - (A*E*t2)/(10*dx) - (6*A*E*w1)/(5*dx^2) - (A*E*t1)/(10*dx), (12*E*I)/dx^3 - (6*A*E*u1)/(5*dx^2) + (6*A*E*u2)/(5*dx^2) + (9*A*E*t1^2)/(70*dx) + (9*A*E*t2^2)/(70*dx) + (108*A*E*w1^2)/(35*dx^3) + (108*A*E*w2^2)/(35*dx^3) + (27*A*E*t1*w1)/(35*dx^2) - (27*A*E*t1*w2)/(35*dx^2) + (27*A*E*t2*w1)/(35*dx^2) - (27*A*E*t2*w2)/(35*dx^2) - (216*A*E*w1*w2)/(35*dx^3),                           (3*A*E*t2^2)/280 - (3*A*E*t1^2)/280 - (6*E*I)/dx^2 + (A*E*u1)/(10*dx) - (A*E*u2)/(10*dx) - (27*A*E*w1^2)/(70*dx^2) - (27*A*E*w2^2)/(70*dx^2) - (3*A*E*t1*t2)/140 - (9*A*E*t2*w1)/(35*dx) + (9*A*E*t2*w2)/(35*dx) + (27*A*E*w1*w2)/(35*dx^2);
+    (A*E*t1)/30 - (2*A*E*t2)/15 - (A*E*w1)/(10*dx) + (A*E*w2)/(10*dx),                                                       (6*E*I)/dx^2 + (3*A*E*t1^2)/280 - (3*A*E*t2^2)/280 - (A*E*u1)/(10*dx) + (A*E*u2)/(10*dx) + (27*A*E*w1^2)/(70*dx^2) + (27*A*E*w2^2)/(70*dx^2) + (3*A*E*t1*t2)/140 + (9*A*E*t2*w1)/(35*dx) - (9*A*E*t2*w2)/(35*dx) - (27*A*E*w1*w2)/(35*dx^2),                                                                        (A*E*u1)/30 - (A*E*u2)/30 + (2*E*I)/dx - (3*A*E*dx*t1^2)/280 - (3*A*E*dx*t2^2)/280 + (3*A*E*t1*w1)/140 - (3*A*E*t1*w2)/140 + (3*A*E*t2*w1)/140 - (3*A*E*t2*w2)/140 + (A*E*dx*t1*t2)/70,             (2*A*E*t2)/15 - (A*E*t1)/30 + (A*E*w1)/(10*dx) - (A*E*w2)/(10*dx),                                                       (3*A*E*t2^2)/280 - (3*A*E*t1^2)/280 - (6*E*I)/dx^2 + (A*E*u1)/(10*dx) - (A*E*u2)/(10*dx) - (27*A*E*w1^2)/(70*dx^2) - (27*A*E*w2^2)/(70*dx^2) - (3*A*E*t1*t2)/140 - (9*A*E*t2*w1)/(35*dx) + (9*A*E*t2*w2)/(35*dx) + (27*A*E*w1*w2)/(35*dx^2), (A*E*dx*t1^2)/140 - (3*A*E*dx*t1*t2)/140 + (3*A*E*t1*w1)/140 - (3*A*E*t1*w2)/140 + (3*A*E*dx*t2^2)/35 - (3*A*E*t2*w1)/140 + (3*A*E*t2*w2)/140 + (9*A*E*w1^2)/(70*dx) - (9*A*E*w1*w2)/(35*dx) + (9*A*E*w2^2)/(70*dx) - (2*A*E*u1)/15 + (2*A*E*u2)/15 + (4*E*I)/dx];
+
+F = [                                                                                                                                                                                                                                                                                                                                                                                       - ((A*E*(18*w1^2 - 36*w1*w2 + 18*w2^2))/30 - (A*E*dx*(30*u1 - 30*u2 - 3*t1*w1 + 3*t1*w2 - 3*t2*w1 + 3*t2*w2))/30)/dx^2 - (A*E*(2*t1^2 - t1*t2 + 2*t2^2))/30;
+    (E*(3360*I*w1 - 3360*I*w2 + 288*A*w1^3 - 288*A*w2^3 - A*dx^3*t1^3 - A*dx^3*t2^3 + 1680*I*dx*t1 + 1680*I*dx*t2 + 864*A*w1*w2^2 - 864*A*w1^2*w2 - 28*A*dx^2*t1*u1 + 28*A*dx^2*t1*u2 - 28*A*dx^2*t2*u1 + 28*A*dx^2*t2*u2 + 108*A*dx*t1*w1^2 + 108*A*dx*t1*w2^2 + 108*A*dx*t2*w1^2 + 108*A*dx*t2*w2^2 + 3*A*dx^3*t1*t2^2 + 3*A*dx^3*t1^2*t2 + 36*A*dx^2*t1^2*w1 - 36*A*dx^2*t1^2*w2 + 36*A*dx^2*t2^2*w1 - 36*A*dx^2*t2^2*w2 - 336*A*dx*u1*w1 + 336*A*dx*u1*w2 + 336*A*dx*u2*w1 - 336*A*dx*u2*w2 - 216*A*dx*t1*w1*w2 - 216*A*dx*t2*w1*w2))/(280*dx^3);
+    (E*(5040*I*w1 - 5040*I*w2 + 108*A*w1^3 - 108*A*w2^3 + 24*A*dx^3*t1^3 - 3*A*dx^3*t2^3 + 3360*I*dx*t1 + 1680*I*dx*t2 + 324*A*w1*w2^2 - 324*A*w1^2*w2 - 112*A*dx^2*t1*u1 + 112*A*dx^2*t1*u2 + 28*A*dx^2*t2*u1 - 28*A*dx^2*t2*u2 + 108*A*dx*t1*w1^2 + 108*A*dx*t1*w2^2 + 6*A*dx^3*t1*t2^2 - 9*A*dx^3*t1^2*t2 - 9*A*dx^2*t1^2*w1 + 9*A*dx^2*t1^2*w2 + 9*A*dx^2*t2^2*w1 - 9*A*dx^2*t2^2*w2 - 84*A*dx*u1*w1 + 84*A*dx*u1*w2 + 84*A*dx*u2*w1 - 84*A*dx*u2*w2 - 216*A*dx*t1*w1*w2 + 18*A*dx^2*t1*t2*w1 - 18*A*dx^2*t1*t2*w2))/(840*dx^2);
+    ((A*E*(18*w1^2 - 36*w1*w2 + 18*w2^2))/30 - (A*E*dx*(30*u1 - 30*u2 - 3*t1*w1 + 3*t1*w2 - 3*t2*w1 + 3*t2*w2))/30)/dx^2 + (A*E*(2*t1^2 - t1*t2 + 2*t2^2))/30;
+    -(E*(3360*I*w1 - 3360*I*w2 + 288*A*w1^3 - 288*A*w2^3 - A*dx^3*t1^3 - A*dx^3*t2^3 + 1680*I*dx*t1 + 1680*I*dx*t2 + 864*A*w1*w2^2 - 864*A*w1^2*w2 - 28*A*dx^2*t1*u1 + 28*A*dx^2*t1*u2 - 28*A*dx^2*t2*u1 + 28*A*dx^2*t2*u2 + 108*A*dx*t1*w1^2 + 108*A*dx*t1*w2^2 + 108*A*dx*t2*w1^2 + 108*A*dx*t2*w2^2 + 3*A*dx^3*t1*t2^2 + 3*A*dx^3*t1^2*t2 + 36*A*dx^2*t1^2*w1 - 36*A*dx^2*t1^2*w2 + 36*A*dx^2*t2^2*w1 - 36*A*dx^2*t2^2*w2 - 336*A*dx*u1*w1 + 336*A*dx*u1*w2 + 336*A*dx*u2*w1 - 336*A*dx*u2*w2 - 216*A*dx*t1*w1*w2 - 216*A*dx*t2*w1*w2))/(280*dx^3);
+    (E*(5040*I*w1 - 5040*I*w2 + 108*A*w1^3 - 108*A*w2^3 - 3*A*dx^3*t1^3 + 24*A*dx^3*t2^3 + 1680*I*dx*t1 + 3360*I*dx*t2 + 324*A*w1*w2^2 - 324*A*w1^2*w2 + 28*A*dx^2*t1*u1 - 28*A*dx^2*t1*u2 - 112*A*dx^2*t2*u1 + 112*A*dx^2*t2*u2 + 108*A*dx*t2*w1^2 + 108*A*dx*t2*w2^2 - 9*A*dx^3*t1*t2^2 + 6*A*dx^3*t1^2*t2 + 9*A*dx^2*t1^2*w1 - 9*A*dx^2*t1^2*w2 - 9*A*dx^2*t2^2*w1 + 9*A*dx^2*t2^2*w2 - 84*A*dx*u1*w1 + 84*A*dx*u1*w2 + 84*A*dx*u2*w1 - 84*A*dx*u2*w2 - 216*A*dx*t2*w1*w2 + 18*A*dx^2*t1*t2*w1 - 18*A*dx^2*t1*t2*w2))/(840*dx^2)];
+
diff --git a/examples/vonKarman/model/Geometry.m b/examples/vonKarman/model/Geometry.m
new file mode 100644
index 0000000..dadfb23
--- /dev/null
+++ b/examples/vonKarman/model/Geometry.m
@@ -0,0 +1,9 @@
+function [Geo] = Geometry(eps,a)
+%% Geometry for a Von Karman Beam element aligned along the x axis (all units SI)
+Geo.eps = eps;
+Geo.l = 1 ; % length 
+Geo.h = eps*Geo.l; % height
+Geo.b = Geo.l/10; % width
+Geo.A = Geo.b*Geo.h; % cross section area
+Geo.I =  Geo.b*Geo.h^3/12; % second moment of cross section
+Geo.a = a; % height of midpoint relative to the ends (measure of curvature of the beam)
diff --git a/examples/vonKarman/model/NonlinearityJacobianVK.m b/examples/vonKarman/model/NonlinearityJacobianVK.m
new file mode 100644
index 0000000..915cea7
--- /dev/null
+++ b/examples/vonKarman/model/NonlinearityJacobianVK.m
@@ -0,0 +1,16 @@
+function [DS] = NonlinearityJacobianVK(model,x)
+% The finite element model is taken from the following article: 
+% Jain, S., Tiso, P., & Haller, G. (2018). Exact nonlinear model reduction 
+% for a von Kármán beam: slow-fast decomposition and spectral submanifolds. 
+% Journal of Sound and Vibration, 423, 195–211. 
+% https://doi.org/10.1016/J.JSV.2018.01.049
+
+q = zeros(model.nDOF,1);
+q(model.freeDOFs) = x;
+
+[Kt,~] = AssembleTangentStiffness(model,q);
+
+DS = Kt - model.K;
+DS = DS(model.freeDOFs,model.freeDOFs);
+end
+
diff --git a/examples/vonKarman/model/NonlinearityVK.m b/examples/vonKarman/model/NonlinearityVK.m
new file mode 100644
index 0000000..166ec62
--- /dev/null
+++ b/examples/vonKarman/model/NonlinearityVK.m
@@ -0,0 +1,14 @@
+function [S] = NonlinearityVK(model,x)
+% The finite element model is taken from the following article: 
+% Jain, S., Tiso, P., & Haller, G. (2018). Exact nonlinear model reduction 
+% for a von Kármán beam: slow-fast decomposition and spectral submanifolds. 
+% Journal of Sound and Vibration, 423, 195–211. 
+% https://doi.org/10.1016/J.JSV.2018.01.049
+
+q = zeros(model.nDOF,1);
+q(model.freeDOFs) = x;
+[~,F] = AssembleTangentStiffness(model,q);
+S = F - model.K * q;
+S = S(model.freeDOFs);
+end
+
diff --git a/examples/vonKarman/model/PlotDeformation.m b/examples/vonKarman/model/PlotDeformation.m
new file mode 100644
index 0000000..ef9fbd2
--- /dev/null
+++ b/examples/vonKarman/model/PlotDeformation.m
@@ -0,0 +1,25 @@
+function PlotDeformation(model,displacement,scale)
+
+l = model.nElements*model.dx;
+u = displacement(1:3:end,:); w = displacement(2:3:end,:);
+dis = sqrt(u.^2 + w.^2);
+maxdispl = max(dis);
+
+sf = scale*l./maxdispl;
+sf = diag(sf);
+if strcmp(scale, 'TRUE')
+    sf = 1;
+end
+% sf = scale;
+
+x = 0:model.dx:l;
+% y = zeros(size(x));
+xd = x' + u*sf;
+yd = w*sf;
+
+% plot(x,y,'-k','LineWidth',2)
+% hold on
+plot(xd,yd)
+% drawnow
+% % axis equal
+% plot(xd(model.nElements/4),yd(model.nElements/04), '*', 'Markersize' , 6)
diff --git a/examples/vonKarman/model/S11.m b/examples/vonKarman/model/S11.m
new file mode 100644
index 0000000..4715b9e
--- /dev/null
+++ b/examples/vonKarman/model/S11.m
@@ -0,0 +1,80 @@
+function X = S11(model)
+% This function assembles the quadratic nonlinearities for the curved beam
+% example.
+
+E = model.E;
+A = model.A;
+dx = model.dx;
+P = cell(1,6); 
+
+P{1} =...
+[ 0,                 0,              0, 0,                 0,              0;...
+  0, -(3*A*E)/(5*dx^2), -(A*E)/(20*dx), 0,  (3*A*E)/(5*dx^2), -(A*E)/(20*dx);...
+  0,    -(A*E)/(20*dx),      -(A*E)/15, 0,     (A*E)/(20*dx),       (A*E)/60;...
+  0,                 0,              0, 0,                 0,              0;...
+  0,  (3*A*E)/(5*dx^2),  (A*E)/(20*dx), 0, -(3*A*E)/(5*dx^2),  (A*E)/(20*dx);...
+  0,    -(A*E)/(20*dx),       (A*E)/60, 0,     (A*E)/(20*dx),      -(A*E)/15];
+
+P{2} =...
+[                 0, -(3*A*E)/(5*dx^2), -(A*E)/(20*dx),                 0,  (3*A*E)/(5*dx^2), -(A*E)/(20*dx);...
+  -(3*A*E)/(5*dx^2),                 0,              0,  (3*A*E)/(5*dx^2),                 0,              0;...
+     -(A*E)/(20*dx),                 0,              0,     (A*E)/(20*dx),                 0,              0;...
+                  0,  (3*A*E)/(5*dx^2),  (A*E)/(20*dx),                 0, -(3*A*E)/(5*dx^2),  (A*E)/(20*dx);...
+   (3*A*E)/(5*dx^2),                 0,              0, -(3*A*E)/(5*dx^2),                 0,              0;...
+     -(A*E)/(20*dx),                 0,              0,     (A*E)/(20*dx),                 0,              0];    
+
+ P{3} = ...
+[              0, -(A*E)/(20*dx), -(A*E)/15,              0,  (A*E)/(20*dx),  (A*E)/60;...
+  -(A*E)/(20*dx),              0,         0,  (A*E)/(20*dx),              0,         0;...
+       -(A*E)/15,              0,         0,       (A*E)/15,              0,         0;...
+               0,  (A*E)/(20*dx),  (A*E)/15,              0, -(A*E)/(20*dx), -(A*E)/60;...
+   (A*E)/(20*dx),              0,         0, -(A*E)/(20*dx),              0,         0;...
+        (A*E)/60,              0,         0,      -(A*E)/60,              0,         0];    
+
+P{4} = ...
+[ 0,                 0,              0, 0,                 0,              0;...
+  0,  (3*A*E)/(5*dx^2),  (A*E)/(20*dx), 0, -(3*A*E)/(5*dx^2),  (A*E)/(20*dx);...
+  0,     (A*E)/(20*dx),       (A*E)/15, 0,    -(A*E)/(20*dx),      -(A*E)/60;...
+  0,                 0,              0, 0,                 0,              0;...
+  0, -(3*A*E)/(5*dx^2), -(A*E)/(20*dx), 0,  (3*A*E)/(5*dx^2), -(A*E)/(20*dx);...
+  0,     (A*E)/(20*dx),      -(A*E)/60, 0,    -(A*E)/(20*dx),       (A*E)/15];  
+
+P{5} = ...
+[                 0,  (3*A*E)/(5*dx^2),  (A*E)/(20*dx),                 0, -(3*A*E)/(5*dx^2),  (A*E)/(20*dx);...
+   (3*A*E)/(5*dx^2),                 0,              0, -(3*A*E)/(5*dx^2),                 0,              0;...
+      (A*E)/(20*dx),                 0,              0,    -(A*E)/(20*dx),                 0,              0;...
+                  0, -(3*A*E)/(5*dx^2), -(A*E)/(20*dx),                 0,  (3*A*E)/(5*dx^2), -(A*E)/(20*dx);...
+  -(3*A*E)/(5*dx^2),                 0,              0,  (3*A*E)/(5*dx^2),                 0,              0;...
+      (A*E)/(20*dx),                 0,              0,    -(A*E)/(20*dx),                 0,              0];
+
+P{6} = ...
+[              0, -(A*E)/(20*dx),  (A*E)/60,              0,  (A*E)/(20*dx), -(A*E)/15;...
+  -(A*E)/(20*dx),              0,         0,  (A*E)/(20*dx),              0,         0;...
+        (A*E)/60,              0,         0,      -(A*E)/60,              0,         0;...
+               0,  (A*E)/(20*dx), -(A*E)/60,              0, -(A*E)/(20*dx),  (A*E)/15;...
+   (A*E)/(20*dx),              0,         0, -(A*E)/(20*dx),              0,         0;...
+       -(A*E)/15,              0,         0,       (A*E)/15,              0,         0];
+
+S = cell(1,model.nDOF);
+for i = 1:model.nDOF    
+    S{i} = zeros(model.nDOF);
+end
+for i = 1:model.nElements 
+    index = 3*(i-1)+1:3*(i+1);
+    Te = compute_rotation_matrix(model.nodes(i:i+1,:));
+    for j = 1:6
+        PPP = zeros(model.nDOF);
+        PPP(index,index) = Te.' * P{j} * Te;
+        ii = index(j);
+        S{ii} = S{ii}+ PPP;
+    end
+end
+
+X = cell(1,length(model.freeDOFs));
+for j = 1:length(model.freeDOFs)
+    i = model.freeDOFs(j);
+    AJ = S{i};
+    X{j} = AJ(model.freeDOFs,model.freeDOFs);
+end
+
+end
\ No newline at end of file
diff --git a/examples/vonKarman/model/compute_rotation_matrix.m b/examples/vonKarman/model/compute_rotation_matrix.m
new file mode 100644
index 0000000..b5b70be
--- /dev/null
+++ b/examples/vonKarman/model/compute_rotation_matrix.m
@@ -0,0 +1,13 @@
+function T = compute_rotation_matrix(nodes)
+T = speye(6,6);
+
+node1=nodes(1,:);
+node2=nodes(2,:);
+l=norm(node2-node1);
+e_x = (node2-node1)/l;
+R = [e_x(1) e_x(2);
+    -e_x(2) e_x(1)];
+
+T(1:2,1:2) = R;
+T(4:5,4:5) = R;
+end
\ No newline at end of file
diff --git a/examples/vonKarman/model/sortVMsVK.m b/examples/vonKarman/model/sortVMsVK.m
new file mode 100644
index 0000000..e36b37c
--- /dev/null
+++ b/examples/vonKarman/model/sortVMsVK.m
@@ -0,0 +1,11 @@
+function [w,VN] = sortVMsVK(V,M,D)
+
+w = (diag(D));
+[w,pos] = sort(w);
+VN = V(:,pos);
+
+
+mu = diag(VN'*M*VN);
+VN = VN*diag(1./sqrt(mu));
+
+end
diff --git a/modeselection/S11reduct.m b/modeselection/S11reduct.m
new file mode 100644
index 0000000..3ccd40f
--- /dev/null
+++ b/modeselection/S11reduct.m
@@ -0,0 +1,16 @@
+function S = S11reduct(SS,X)
+% This function converts the quadratic nonlinearity coefficients from
+% their general form to the desired form in modal coordinates
+    U = SS.U;
+    U2 = SS.U2;
+    B = size(U2);
+    S = cell(1,B(2));
+    n = B(2)+SS.n;
+    
+    for j = 1:B(2)
+        S{j} = zeros(SS.n,SS.n);
+        for i = 1:n
+            S{j} = S{j} + U2(i,j) * U.' * X{i} * U;
+        end
+    end    
+end
\ No newline at end of file
diff --git a/modeselection/copyfromfull.m b/modeselection/copyfromfull.m
new file mode 100644
index 0000000..22b2f69
--- /dev/null
+++ b/modeselection/copyfromfull.m
@@ -0,0 +1,15 @@
+function class = copyfromfull(class,full)
+% copies properties of the second entry to the first entry
+if ~isempty(full.f)
+    class.f = full.f;
+end
+if ~isempty(full.T)
+    class.T = full.T;
+    class.n_steps = full.n_steps;
+end
+if ~isempty(full.order)
+    class.order = full.order;
+end
+    class.S = full.S;
+    class.DS = full.DS;
+end
\ No newline at end of file
diff --git a/modeselection/modeselect.m b/modeselection/modeselect.m
new file mode 100644
index 0000000..47ce78b
--- /dev/null
+++ b/modeselection/modeselect.m
@@ -0,0 +1,101 @@
+function modes = modeselect(full,param,X) 
+%% Automated nonlinear mode selection procedure: described in Algorithm 1 of
+% the following article 
+% G. Buza, S. Jain, G. Haller, Using Spectral Submanifolds for Optimal 
+% Mode Selection in Model Reduction, (2020) Preprint available on arXiv.org
+
+
+try
+    modes = param.initialmodes; 
+catch    
+    % initial set based on spectral quotient analysis
+    realpart = min(abs(real(-full.zeta.*full.omega0 + sqrt(full.zeta.^2-1).*full.omega0)),...
+        abs(real(-full.zeta.*full.omega0 - sqrt(full.zeta.^2-1).*full.omega0)));
+    [rpsorted,ind] = sort(realpart);
+    ratio = zeros(1,length(rpsorted)-1);
+    for i = 1:length(rpsorted)-1
+        ratio(i) = rpsorted(i+1)/rpsorted(i);
+    end
+    cut = param.nmodes;
+    while cut > ceil(param.nmodes/3)    % limit the number of modes selected via s.q.a.
+        [~,cut] = max(ratio);
+        if cut ~= 1
+        ratio = ratio(1:cut-1);
+        end
+    end
+    modes = ind(1:cut).';
+
+    % refinement based on linear response to forcing (if given)
+    if ~isempty(full.f)
+        [x_linf,~] = full.LinearResponse();
+        z = full.U.' * full.M * x_linf;
+        b = size(z);
+        NORM = zeros(1,b(1));
+        for j = 1:b(1)
+            NORM(j) = norm(z(j,:));
+        end
+        while length(modes) < ceil(param.nmodes*2/3)    % limit the number of modes selected via linear response   
+            if sum(NORM(modes)) > 0.9*sum(NORM)
+                break
+            end
+            modesC = setdiff(1:full.n,modes);
+            [~,ind2] = max(NORM(modesC));
+            modes = [modes ind2];
+        end
+        percent = 100*sum(NORM(modes))/sum(NORM);
+        disp(['The initial mode selection recovers ' num2str(percent) '% of the linear response.'])
+    else
+        disp('Since no forcing was provided for the SSR class, it is not guaranteed that the mode selection recovers the linear response.')
+    end
+end
+disp('Initial linear selection:')
+disp(modes)
+
+% nonlinear selection
+itnumber = 0;
+while length(modes) < param.nmodes
+    red = SSR(full.M,full.C,full.K,modes);
+    red = copyfromfull(red,full);       % copies general props
+    switch param.example                % obtaining the nonlinear coefficients is problem dependent
+        case 'simple'                   % heavy on memory
+            R = S11reduct(red,X);
+        case 'plate'                    % memory preserving approach
+            V = zeros(param.fulldof,red.n);
+            VC = zeros(param.fulldof,length(red.M)-red.n);
+            V(param.freedofs,:) = red.U;
+            VC(param.freedofs,:) = red.U2;
+            [K2,~] = AssembleTensors(param.model,V,VC); % computes quadratic coeffs R_i
+            R = cell(1,size(red.U2,2));
+            for i = 1:size(red.U2,2)
+                R{i} = K2(i,:,:);
+            end
+    end
+    W = red.SSM2(R);                    % computing SSM
+    A = [zeros(red.n),eye(red.n);-red.K1,-red.C1];
+    [~ ,scaldir] = scalg(A,W);          % computing directional curvatures
+    [scnorm,scindex] = sort(abs(scaldir));
+    Jmodes = setdiff(1:full.n,modes);
+    scindex = fliplr(scindex);
+    crit = sum(scnorm) * (1-param.curvtolerance);
+    Z = 0;
+    numrec = 0;
+    while Z < crit 
+        Z = Z + max(scnorm);
+        scnorm = setdiff(scnorm,max(scnorm));
+        numrec = numrec + 1;
+    end
+    recomm = Jmodes(scindex(1:numrec)); % recommnended modes
+    check = param.nmodes - (length(modes) + length(recomm)); % check if it overshoots the number of allowed modes
+    if check < 0
+        modes = sort([modes recomm(1:(param.nmodes - length(modes)))]);
+    else
+        modes = sort([modes recomm]);
+    end
+    disp(modes)
+    itnumber = itnumber + 1;
+    if ~param.type
+        break
+    end
+end
+disp(['The nonlinear mode selection procedure was performed ' num2str(itnumber) ' time(s).'])
+end
\ No newline at end of file
diff --git a/modeselection/scalg.m b/modeselection/scalg.m
new file mode 100644
index 0000000..919105c
--- /dev/null
+++ b/modeselection/scalg.m
@@ -0,0 +1,24 @@
+function [scal,scaldir] = scalg(A,W)
+% This function computes directional scalar curvatures.
+
+m = length(W{1});
+
+SUM = 0;
+scaldir = zeros(1,length(W));
+for l = 1:length(W)
+    Wl = W{l};
+    sumk = 0;
+    for i = 1:m
+        for j = 1:m
+            a = 4 * ( Wl(i,i) * Wl(j,j) - (Wl(i,j))^2 +...
+                4* (Wl(i,:)*A(:,i))*(Wl(j,:)*A(:,j))-...
+                (Wl(i,:)*A(:,j) + Wl(j,:)*A(:,i))^2);
+            SUM = SUM + a;
+            sumk = sumk+a;
+        end
+    end
+    scaldir(l) = sumk;
+end
+
+scal = SUM/(m^2);
+end
\ No newline at end of file