my matlab code is this for 8.pdf file:

%Design and Optimization of a MEMS Electret-Based Capacitive Energy Scavenger

%coding date start;16/05/1403

%majid akbari

%--------------------------------------------------------------------------

clc;

clear all;

Vel = 150; % surface potential voltage of electret layer [V]

Lx = 5e-3; % length of device in x direction [m]

Ly = 2.3e-3; % length of device in y direction [m]

b = 100e-6; % height of fingers [m]

B = 25.5e-6; % damping coeffient [Ns/m]

Cel = 16.8e-12; % capacitance of electret layer [F]

freq = 810; % freguency

amp = 5e-6; % excitation amplitiude [m]

h = 140e-6; % initial overlap between fingers [m]

rho = 2330; % material density of si [Kgr/m^3]

epsz = 8.854e-12; %vaccum permitivity [F/m]

R = 70e+6; % electrical load

%--------------------------------------------------------------------------

ar = 25; % aspect ratio of system

delta = b/ar; % the gap between movable and fixed fingers [m]

d = b/ar; % width of fingers [m]

H = h * 2; % length of fingers [m]

w = h * 2; % width of fixed electrode [m]

omega1 = 2*pi*freq; % angular freq [rad/s]

A = amp * (omega1^2); % accelaration of excitation source [m/s^2]

W = Lx - 2 * (2 * H + w - h); % width of mass [m]

N = ceil(Ly / (2 * d + 2 * delta)); % number of fingers

eps = 2;

Cz = 2 * N * epsz * b * h / delta; % initial capacitor

Vz = Cel * Vel / (2 * Cz + Cel); % quiescent voltage for each variable capacitor

mass = rho * b * Ly * (Lx- 9 * h); % mass [kg]

%--------------------------------------------------------------------------

%normalization parameter

Czero = (epsz * Ly * b * h) / (2 * d^2);

ETA = (Cz * Vz^2) / (mass * h^2 * omega1^2);

ZETA = B / (2 * mass * omega1);

a = @(t) -A * sin(omega1*t); %acceleration caused on the device by the external vibration source

ap = @(t) a(t) / (h * omega1^2); %a^(t)

omega2p = 1 / (R * Cz * omega1);

omegaelp = 1 / (R * Cel * omega1); %omega electret wel^(t)

Ielp = Vel / (R * Cz * Vz * omega1);%Iel^

Qzero = Cz * Vz;

k = mass * (omega1)^2;

%--------------------------------------------------------------------------

%odefun= @(t,y) [-y(1)*delta / (2 * N * R * epsz * b *(h-y(3))) - (y(1) + y(2)) / (Cel*R) + Vel/R;

% -y(2)*delta / (2 * N * R * epsz * b *(h+y(3))) - (y(1) + y(2)) / (Cel*R) + Vel/R;

% y(4);

% (y(2)^2*delta) / (4*N*mass*epsz*b*(h+y(3))^2)-(y(1)^2*delta)/(4*N*mass*epsz*b*(h-y(3))^2) - (k/mass) * y(3) - (B/mass) * y(4) -a(t)];

%[t,y]=ode15i(odefun,[0 30],[0.82*Qzero 0.9*Qzero delta 0]);

%--------------------------------------------------------------------------

%ODE45 FUNCTION

odefun= @(t,y) [-omega2p * (y(1)) / (1 - y(3)) - omegaelp * (y(1) + y(2)) + Ielp;

-omega2p * (y(2)) / (1 + y(3)) - omegaelp * (y(1) + y(2)) + Ielp;

y(4);

-(ETA / 2) * (y(1)/(1 - y(3)))^2 + (ETA / 2) * (y(2)/(1 + y(3)))^2 - (y(3)) - (2 * ZETA * y(4)) - ap(t)];

opts = odeset('RelTol',1.323489e-13,'AbsTol',1e-27);

[t,y] = ode45(odefun,[0,5],[0.82 0.9 -0.7 0],opts);

%-------------------------------------------------------------------------

plot(t, y(:,1))

hold on;

plot(t, y(:,2),"green")

hold on;

plot(t, y(:,3),"red")

but i get this result in output in file 11.png