nonlinear ode to solve RMS power gives me fluctuating trend and zero values!!

Hi,

I have a non-linear mechanical system, I am solving the power generated from the movement of the vehicle suspension system

I used simulink to solve the non-linear ode and extract the relative velocity to use it on the power equation.

I think my code is correct but don't know why I am getting zero values for the RMS power at some points. I should be getting

increasing trendline on my results figure! I noticed that changing the time steps play significate role in chaning the results, similarly when I change the variance for the white noise.

my code works fine when i solve for other parameters, but doesn't work good when i solve with respect to velocity changing overtime!

is there any why to make the graph looks more smoother and the code run faster?

clc
clear
close all


%% ((@@ calculating for power w/ respect to space Velocity changing @@)) 

%   =================================================================


%% /// parameters for Piezoelectric harvester

T= 30;                                % period                                         (Sec)
l=0.1;                                % length of the magnetic anf Piezo                (m)
ld=0.01;                              % lead of the ball screw                          (m)
w=0.002;                              % wedith of the Piezo & magnetic                  (m)
d=0.0005;                             % spacing btw the rotator & stator                (m)
r1=0.5;                               % inner raduis of the stator ring                 (m)
r2=0.05;                              % outer raduis of the rotator ring                (m)
n1= 0;                                % initial value for suspension velocity           (m)
Br= 1.5;                              % residual flux density of magnetic               (T)
n2=abs(pi*r2/w);                      % #of magnet slabs only on rotator ring
d33= 6.4e-10;                         % piezoelectric coefficient                      (CN^-1)
num_trying=12;                        % final number of RMS-power we need to calculate
tm = 0.01;                            % thickness of magnetic slab                      (m)
tp = 0.01;                            % thickness of Piezo patches                      (m)
d0=0.001;                             %                                                 (m)
Cv_dash=0.375*10^-9 ;                 % Cv unit changed to farad                       (Farad)
Cv=Cv_dash*l*w*0.0001/(0.01*0.01*tp); % electric capacity of piezoe patch              (Farad)


%% /// parameters for vehicle and road 

mi=5.3;                               % is inertial mass                    (kg)
mb=362;                               % is the body mass                    (kg)
mt=30;                                % is the wheel mass                   (kg)
c=1.4*1000;                           % damper                              (N.s/m)
k=20.1*1000;                          % spring                              (N/m)
kt=182.1*1000;                        % Tire stiffness                      (N/m)
Gq=1024e-6;                           % road roughness coefficient           m^3
no=0.1;                               % reference spatial frequency          m^-1
n00=0.011;                            % minimal boundary frequency           m^-1                    

v=linspace((10*1000)/3600,(120*1000)/3600,num_trying); % speed of the vehicle        m/s

V=linspace(10,120,num_trying);


%% Numerical analysis 

Pe= 0;
Pe_rms_arr=zeros(num_trying,1);


for s=1:1:num_trying

Bd=(Br/pi)*((atand((l*w)/(2*d*sqrt(4*(d)^2+(l)^2+(w)^2))))-(atand((l*w)/(2*(tm+d)*sqrt(4*(tm+d)^2+(l)^2+(w)^2)))));
fd=(1.749+1.145*exp(-d/d0))*(10^6);
Fm=l*w*((tm)^1/3)*Br*abs(Bd)*fd;
Fp=2*((d33)^2*(Fm)^2*(n2)^2)/Cv*ld;   



%% /// (((( Run the Simulink model using the 'sim' command )))) 
sim('nonlinear_model_v.slx')

%% /// extracting Data generated by Simulink

 ZdotData = ans.z_dot;
 Zdot     = ZdotData.Data;    % Z'= (Zb' - Zt')
 n1       = Zdot/ld;          % n1= (zb' - zt') / ld 
 
%% Calculating the power generation

y= 1;
time_step=100*T;
delta_t=T/time_step;
Pe_arr=zeros(time_step,1);

for t=0:delta_t:T
     
    %   loop for all peizoelectric patches  for one time step
    
       for N2=1:1:2*n2
           for s2=1:1:length(Zdot)
          % Pe=Pe+(((d33)^2*n1(s2)*N2*pi*(Fm)^2*abs(sind(n1(s2)*N2*pi*t)*cosd(n1(s2)*N2*pi*t)))/Cv) ;
           Pe=Pe+(((d33)^2*((Zdot(s2)/ld)*N2*pi*(Fm)^2*abs(sind((Zdot(s2))/ld)*N2*pi*t)*cosd(((Zdot(s2))/ld)*N2*pi*t)))/Cv) ;
           end
           
       end
       
    Pe_arr(y)=Pe; %  store  the power here for each time step
    Pe=0;
    y=y+1;
   
  end  
   %  now calculate the RMS power we got
  %   loop for all power  for j step 
  for j=2:1:time_step
%     *************  sum=Pe_arr(j)^2+Pe_arr(j-1)^2;
      sum=Pe_arr(j)^2-Pe_arr(j-1)^2;
  end
  Pe_rms=sqrt((delta_t/(2*(T-delta_t)))*sum);

  Pe_rms_arr(s)=Pe_rms;
%asd=s

end 



%% figures 

figure 
plot(ans.z_dot)
xlabel('Time/s')
ylabel('zdot(t)') 
grid on

figure 
plot(ans.road)
xlabel('Time/s')
ylabel('q(t)') 
grid on

figure

plot(V,real(Pe_rms_arr),'b--*')
title('the RMS power versus with driving speed');
ylabel('RMS of the power (Watt)');
xlabel('speed (Km/h)');
legend('The RMS')
grid on