ODE45 is not solving (not getting triggered) the differential equation

Hello everyone,

I am using ODE45 to solve a large number of second order differential equation (Right hand side of this ODE is non zero function (a_i(i))) in my code. In my code ODE45 is not solving the differential equation. It is generating the output vector same as initial condition (displacement). I am getting Output (displacement) at all the time step same as initial displacement (i.e., it is directly giving initial displacement as final output for all ODE, which is not correcct solution). But output should be a sinusoidal function. I have posted my code below. Please please please help me to solve this problem.

Thanks in advance.

%% AMM mode 1 ODE dsolve
%% ================================================================
clc; clear all; close all
syms eta1(t)
m1=1; c1=0; k1=9.86;
Deta1 = diff(eta1);
ode = m1*diff(eta1,t,2) + c1*Deta1 + k1*eta1 == 0 ;
cond1 = eta1(0)== 8/(pi)^3; cond2 = Deta1(0)== 0;
conds = [cond1 cond2];
eta1(t) = dsolve(ode,conds)
syms etaa1
etaa1 = eta1(t);               % eta1(t)
%% AMM mode 3 ODE dsolve
%% ============================================================
clc; clearvars -except keepVariables etaa1; 
close all
syms eta3(t)
m3=1; c3=0; k3=88.82;
Deta3 = diff(eta3);
ode = m3*diff(eta3,t,2) + c3*Deta3 + k3*eta3 == 0; 
cond1 = eta3(0)== 8/(3*pi)^3; cond2 = Deta3(0)== 0;
conds = [cond1 cond2];
eta3(t) = dsolve(ode,conds)
syms etaa3
etaa3 = eta3(t);               % eta3(t)
%% AMM mode 5 ODE dsolve    
%% =======================================================
clc; clearvars -except keepVariables etaa1 etaa3;
close all
syms eta5(t)
m5=1; c5=0; k5=246.74;
Deta5 = diff(eta5);
ode = m5*diff(eta5,t,2) + c5*Deta5 + k5*eta5 == 0; 
cond1 = eta5(0)== 8/(5*pi)^3; cond2 = Deta5(0)== 0;
conds = [cond1 cond2];
eta5(t) = dsolve(ode,conds);
syms etaa5
etaa5 = eta5(t);                  % eta5(t) 
%% displacement expression due to three modes
syms x
u(x,t) = etaa1*sin(pi*x)+etaa3*sin(3*pi*x)+etaa5*sin(5*pi*x);

%%
velo(x,t) = diff(u,t);

%% divide length elmentwise
 x_i = 0:0.1:1;                      % Locaion of points 
 for i = 1:length(x_i)
    u_i(i) = subs(u,x,x_i(i));       % Displacement of points
      
 end
u_i;                                 % Displacement of points
v_i = diff(u_i,t);                   % Velocity of points
a_i = (diff(diff(u_i,t)));           % Acceleration of points

%% Initial condition vectors for all elements
IDvec = subs(-u_i,t,0);    % initial  displacement vector for points
IVvec = subs(v_i,t,0);    % initial  velocity vector for points 
%%    OK upto this point 
global k Th0 m c 
e=1;  m=1;  c=0;  k=0; 

t_vector = [] 
ui_vector = []
   
no_of_elements = length(x_i);
for i = 1:no_of_elements   % element for which ODE is being solved for contact surface
%==============================================
Thid = double(IDvec(i))	% position from where the mass is released
Thiv = 0	% inital velocity (at t=0)

Th0 = 0.005	% Position of the wall 

tfi = 5	% time final time    

% ode45 (function call)====================================
ui0=[Thid Thiv]     % Initial conditions
tspan=[0 tfi]
options = odeset('Events',@aakash_funct); % event 

t_vec = []
ui_vec = []

tcon=0      

 ii=1;
while tcon<tfi		%% to check whether final time is reached or not
        tspan = tcon:0.01:(tfi+0.01);
       [t,ui] = ode45(@mode123dsolve,tspan,ui0,options); % ode45 solver=========

    
    ui0=[ui(end,1)  -(e*ui(end,2))]	%% modified initial conditions after impact
    tcon=t(end,1)
   
     
   t_vec = [t_vec;t]
   ui_vec = [ui_vec;ui] 
end
%     t_vec = [t_vec;t] 
%     ui_vec = [ui_vec;ui]

   t_vector = [t_vector t_vec] 
   ui_vector = [ui_vector ui_vec] 
end 
 

The function file is given here..

%% Function FILE CODE

function dui=mode123dsolve(t,ui)

global  k Th0 m c a_i i
%% Equations of motion in state-space form
E=[0 1;-k/m -c/m];
F=[0; a_i(i)];        % a_i(i) is a function of time that will be obtained from the script file. 
dui=E*ui+F;

event file code is as follows..

%% Event condition code

function [value,isterminal_1,direction_1] = aakash_funct(t,ui)
global  L Th0  e

value = ui(1)-Th0 ;  % when value = 0, an event is triggered
isterminal_1 = 1;  % terminate after the first event
direction_1 = +1 ;  % get all the zeros