How to set time in matlab coding

Hello everyone,

Hello everyone,

Iam designing piezoelectric energy harvester and Iam using analytical modelling technique for determining voltage using matlab coding. While trying to execute it, iam getting voltage output, but the oscillations are not proper, and the time is increasing randomly upto 8000 seconds. I need to set the voltage [y axis] within -8 millivolts to +8 millivolts and time[x axis] within 0 to 10 milliseconds. I couldn't set the voltage and time in coding. I have attached my coding here. Please, tell me where i need to change in the coding to get correct output. Thank you..

clear all;
close all;

%Beam properties

rho=2700;                  %beam density
Ebeam=70e9;                %young's modulus of beam
Lbeam=0.003;               %beam length
width=0.001;               %beam and piezo width
Lf=0.200;                  %length to forcing f(t)
thickbeam=0.0002;          %beam thick
Area=thickbeam*width;      %area
I=width*thickbeam^3/12;    %Inertia

%Piezo peoperties
g31=-9.5e-3;               %piezoelectric coefficient
Epiezo=50e9;               %young's modulus of piezo
thickpiezo=0.00005e-3;     %piezo thick
Lpiezo=0.003;              %piezo length
d31=-320e-12;              %m/volt dielectric constant
Rsource = 330;          %ohm source resistance
c=(Ebeam*I/(rho*Area))^.5;

%Natural Frequencies

B = 0:.01:50; %(wn/c)^2;
cheq = cos(Lbeam*B).*cosh(Lbeam*B);
plot(B,cheq)
axis([0 30 -5 5])
hold on;
ezplot('-1',[0 30 -5 5])
axis([0 30 -5 5])
xlabel('Beta values');title('Graph of Characteristic equation')
beta = [3.3604;8.4123;14.0766;19.71;25.335];
wn=beta.^2*c;
fn=wn/(2*pi);

%ModeShapes

syms x
for i=1:5,
   ModeShape(i) = cosh(beta(i)*x)-cos(beta(i)*x)-(sinh(beta(i)*Lbeam)...
-sin(beta(i)*Lbeam))/(cosh(beta(i)*Lbeam)+cos(beta(i)*Lbeam))...
*(sinh(beta(i)*x)-sin(beta(i)*x));
end

%Harmonic Forcing Function

w=wn;
zeta=0.03;
wd=wn*(1-zeta^2)^.5;
syms t tt
Force=(0.2)*sin(w(1)*tt)*1/(rho*Area)*subs(ModeShape,x,Lf);
for i=1:5
    q(i)=1/wd(i)*exp(-zeta*wn(i)*t)...
        *int(Force(i)*exp(zeta*wn(i)*tt)*sin(wd(i)*(t-tt)),tt,0,t);
end

%Develop voltage equations

for i=1:5
    before_w3(i)=q(i)*ModeShape(i);
end

w3=sum(before_w3);
K_x_t=vpa(diff(w3,x,2),5);
Kavg=vpa((1/Lpiezo)*int(K_x_t,x,0,Lpiezo),5);
tau = 3:0.01:10;
Moment=Ebeam*I.*Kavg;
psi=(Ebeam*thickbeam)/(Epiezo*thickpiezo);
T=thickbeam/thickpiezo;
Veb_sym=(6*g31*psi*(1+T))/(width*thickpiezo*(1+psi^2*T^2+2*psi*(2+3*T+2*T^2))).*Moment;
Veb=subs(Veb_sym,'t',tau);

plot(Veb,'g');
xlabel('time');
ylabel('voltage');