How to determine convergence of ode45

Hi all,

I was wondering how it is possible to determine convergence of ode45 for a simple damped mass-spring system. So here is the code for the system:

clear
clc
% Inputs
    % Degrees of Freedom
     N = 3 ;

    % Masses (kg)
     m = repmat(100, 1, N) ;

    % Spring coefficients (N/m)
     k = repmat(15, 1, N+1) ;
    
    % Damping coefficients (Ns/m)
    c = repmat(5, 1, N+1) ;
    
    % Maximum force (N)
    f0 = repmat(100, 1, N) ;
    
    % Phase angle between force (rad)
    phi = linspace(0, pi, N) ;
    
    % Frequency of force (Hz)
    freq = 0.2 ;
    
    % Time of simulation (s)
    tmax = 100 ;
    
    % Mass initial positions (m)
    xi = zeros(1, N) ;

% ODE Inputs
    % Degrees of Freedom
    N = length(m) ;
    
    % Mass Matrix
    M = diag(m) ;
    
    % Stiffness Matrix
    k1 = diag(k(1:end-1) + k(2:end)) ;
    k2 = diag( -k(2:end-1), 1) ;
    k3 = diag( -k(2:end-1), -1) ;
    K = k1 + k2 + k3 ;
    
    % Damping Matrix
    c1 = diag(c(1:end-1) + c(2:end)) ;
    c2 = diag( -c(2:end-1), 1) ;
    c3 = diag( -c(2:end-1), -1) ;
    C = c1 + c2 + c3 ;
        
% ODE Initial Conditions
    % Time
    tspan = [0 tmax] ;
    
    % Displacement and velocity
    x0 = [ xi zeros(1, N) ] ;
    
    % ODE Options
    options_set = odeset('Mass', [eye(N) zeros(N) ; zeros(N) M]) ;

    % Solution
    [t, xSol] = ode45(@(t, x) odefcn(t, x, K, C, f0, freq, phi, N), tspan, x0, options_set) ;
 
    % Results
    x = xSol(:, 1:N) ;          % Displacement (m)
    xdot = xSol(:, N+1:end) ;   % Velocity (m/s)

% ODE Function
function dxdt = odefcn(t, x, K, C, f0, freq, phi, N)

% Indices
m = N + 1 ;
n = N * 2 ;

% Preallocate
dxdt = zeros(n, 1) ;

% ODE Equation
dxdt(1:N) = x(m:n) ;
dxdt(m:n) = - K * x(1:N) - C * x(m:n) + f0' .* sin(2*pi*freq * t + phi') ;

end

Now, as I understand, in order to determine convergence I will have to use an EventsFcn. However what perplexes me is how to express the convergence formula within the odefcn, as I can't see how I can compare the one step with the next.

Kind Regards,

KMT.

[Aquatris]

Since it is a simple mass-spring-damper system, why don't you just check the pole locations (eigenvalues of A matrix, A =[zeros(n) eye(n); -M\K -D\K]). If all poles are on the left half plane (negative real parts) than the ODE will converge.

[KostasK]

Yes, this specific problem is simple enough such that this approach could suffice. I think for my above question, its better to calculate the Residuals than the error convergence (which should happen anyway since the solver does not return any warnings/errors).