There is a problem in creating the legend string. I have modififed the code, now error is rectified. But your code is messy. You have to make many changes.
%% ------Energy Harvesting using Piezo Material with Cantilever Structures------
% This promgram intend to use the Continuous Vibration Equation for estimate
% the power generation of a cantilever fully covered by piezoelectric material
% The Math background is provided by the paper "A Distributed Parameter
% Electromechanical Model For Cantilevered Piezoelectric Energy Harvesters" by
% A.Erturk and D.J.Inman.
clear ; close all; clc
%----------------------------INPUT PARAMETERS----------------------------------%
% This simulation can range diffente input parameters and then plot it. Choose
% one parameter and replace by MRANGE(c,2) variable. The MRANGE(c,2) will change
% According to the RANGE Variable bellow
% RANGE MODULE
RANGE = 0.1:0.1:0.5; % Range for used for change the parameters (MAX: 8)
VRANGE = 1:length(RANGE); % Vector with same size of the RANGE;
MRANGE = [VRANGE',RANGE']; % Matrix with The the Index of RANGE
F = 1:0.1:100; % Frequency Range (Hz)
RANGEV0 = ones(length(F),length(RANGE)); % Matrix for Voltage simulation data
RANGEI0 = ones(length(F),length(RANGE)); % Matrix for Current simulation data
RANGEP0 = ones(length(F),length(RANGE)); % Matrix for Power simulation data
% MECHANICAL PARAMETER INPUT
% Beam Support
for c = VRANGE
leB = MRANGE(c,2); % Lenght of the Beam (m)
wiB = 0.02; % Widht of the Beam (m)
thB = 0.0005; % Thickness of the Beam (m)
yoB = 100000000000; % Young's Modulos of the Beam kgf/m2
deB = 7165; % Material Density (kg/m3)
% Piezoelectric material
leP = leB; % Lenght of the Piezo (m) OBS: SHOULD BE SAME AS THE BEAM
wiP = wiB; % Widht of the Piezo (m)
thP = 0.0004; % Thickness of the Piezo (m)
deP = 7800; % Material Density of Piezo (kg/m3)
yoP = 66000000000; % Young's Modulos of the Piezo
% ELECTRICAL PARAMETER INPUT
% Piezoelectric material
d31 = -190*(10^-12); % Piezoelectric Constant at 31 Direction
e33 = 15.93*(10^-9); % Piezoelectric Dielectric Constant (F/m)
% BASE HARMONIC VIBRATION
Y0 = 0.000001; % Peak Displacement of the base (m)
f = F; % Frequency Range (Hz)
w = 2*pi()*f; % Angular frequency (Hz)
end
% ---------------------------FIRST CALCULATIONS--------------------------------%
% CALCULATING THE MASS DENSINTY
% The mass density is the total mass of the system dived by the total length of
% the Beam
volB = thB*wiB*leB;
volP = thP*wiP*leP;
maB = deB*volB;
maP = deP*volP;
massden = (maB+maP)/leB; % (kg/m)
% LOAD RESISTOR
% Resistor conected to the Piezo
rl = 1000000 ;% Ohms
% POSITION RELATIVE TO NEUTRAL AXIS
% The procedure of finding the position of the neutral axis of a composite cross
% section is transform the cross section to a homogeneous cross section of
% single Young’s modulus (Described in Timoshenko and Young Reference)
nY = yoB/yoP; % Ratio of Young's Modulos
hpa = ((thP ^ 2) + (2*nY*thP*thB) + (nY*(thB ^ 2)))/(2*(thP + (nY*thB))); % Distance from the top of the Piezo Layer to the Neutral Axis
hsa = ((thP ^ 2) + (2*thP*thB) + (nY*(thB ^ 2)))/(2*(thP + (nY*thB))); % Distance from the Bottom of the Piezo Layer to the Neutral Axis
hpc = (nY*thB*(thP + thB))/(2*(thP + (nY*thB))); % Distance from the Center of the Piezo Layer to the Neutral Axis
ha = -hsa; % Position of the bottom of the substructure layer from the neutral axis
hb = hpa - thP; % Position of the bottom of the Piezo layer from the neutral axis
hc = hpa; % Position of the top of the layer from the neutral axis
% BENDING STIFFNESS
yi = wiB*(((yoB*((hb^3)-(ha^3))) + (yoP*((hc^3) - (hb^3))))/3);
% COUPLING TERM
coup = - (((yoP*d31*wiB) / ((2*thP))*((hc^2) -(hb^2))));
% UNDAMPED NATURAL FREQUENCY OF THE BEAM
% The sigma constant for Cantilever Structuree is provide by the literature.
% The Solution for sigma is solving the equation 1 + cos (sig) cosh (sig) = 0
% For the three first modes we have:
sig1 = 1.875104069;
sig2 = 4.694091133;
sig3 = 7.854757438;
sig = [sig1;sig2;sig3];
wr1 = (sig1^2) * ((yi / (massden*(leB^4)))^0.5);
wr2 = (sig2^2) * ((yi / (massden*(leB^4)))^0.5);
wr3 = (sig3^2) * ((yi / (massden*(leB^4)))^0.5);
wr = [wr1;wr2;wr3];
% EIGENFUNCTIONS UNDAMPED FREE VIBRATION
ho = ((cosh(sig) + cos(sig))./(sin(sig) + sinh(sig)));
x = leB;
var = ((sig*x)./leB);
eigeinf = ((1/(massden*leB))^0.5).*(cosh(var) - cos(var) - (ho.*(sinh(var)-sin(var))));
% BENDING SLOPE (DERIVATE OF EIGENFUNCTIONS)
bendS = ((1/(massden*leB))^0.5)*(((sig./leB).* sinh(var)) + ((sig./leB).* sin(var)) - (ho.*(((sig./leB).*cosh(var))-((sig./leB).* cos(var)))));
% INTEGRATION OF THE EIGENFUNCTIONS
inteia = ((1/(massden*leB))^0.5)*((sinh(var)./(sig./leB)) - (sin(var)./(sig./leB)) - (ho.*((cosh(var)./(sig./leB))+(cos(var)./(sig./leB)))));
intei0 = ((1/(massden*leB))^0.5)*((sinh(0)./(sig./leB)) - (sin(0)./(sig./leB)) - (ho.*((cosh(0)./(sig./leB))+(cos(0)./(sig./leB)))));
intei = inteia - intei0;
% TIME CONSTANT OF THE PIEZO REPRESENTATION CIRCUIT
tc = (rl*e33*wiP*leP) / thP;
% MODAL CONSTANT TERM
mco = -(((d31*yoP*hpc*thP)/(e33*leB)) * bendS);
% MODAL COUPLING TERM
mc = coup * bendS;
% MODAL DAMPING (ASSUMING DAMPING RATIOS ZETA 1 = 0.01 AND ZETA 2 = 0.013)
ca = 0.654961133;
csI = 1.216501665*(10^-6);
zeta = ((csI*wr)/(2*yi))+ (ca./(2*wr*massden));
% VOLTAGE ESTIMATION
modes = 3;
fs = size(f);
numV = ones(modes,fs(2));
denV = ones(modes,fs(2));
for b = 1:modes
numV(b,:) = ((-1i*massden*w*mco(b)*intei(b))./(((wr(b)^2)-(w.^2)) + (1i*2*zeta(b)*wr(b)*w)));
denV(b,:) = ((1i*w*mc(b)*mco(b))./(((wr(b)^2)-(w.^2)) + (1i*2*zeta(b)*wr(b)*w)));
end
num = numV(1,:) + numV(2,:) + numV(3,:);
den = denV(1,:) + denV(2,:) + denV(3,:);
denV2 =((1 + (1i*w*tc))/tc);
v = num./(den+denV2);
V0 = abs(v);
V = V0.*((w.^2)*Y0);
RANGEV0(:,MRANGE(c,1)) = V0';
RANGEV(:,MRANGE(c,1)) = V';
% CURRENT ESTIMATION
j = (num./((den+denV2)*rl));
I0 = abs(j);
I = I0.*((w.^2)*Y0);
RANGEI0(:,MRANGE(c,1)) = I0';
RANGEI(:,MRANGE(c,1)) = I';
% POWER ESTIMATION
p = v.*j;
P0 = abs(p);
P = P0.*(((w.^2)*Y0).^2);
RANGEP0(:,MRANGE(c,1)) = P0';
RANGEP(:,MRANGE(c,1)) = P';
% PLOTTING DATA
% Input the legend, title, and lables for Axis. Assign the Range variable
% VOLTAGE, CURRENT or POWER for plot (Use RANGEV0, RANGEI0 or RANGEP0 respectively).
xlab = "Frequency (Hz)";
ylab = "FRF Voltage (V / Base Acceleration) ";
tit = "FRF Voltage Variation with the Length(m)";
legS = "Length ";
for Range = RANGEV0
figure(1);
colororder = get (gca, "colororder");
end
labels = cell(size(VRANGE)) ;
for c = VRANGE
hold on
h = plot (f,Range(c,:));
hold on;
set (h, "color", colororder(c,:));
labels{c} =strcat('Length = ',num2str(MRANGE(c,2)));
hold off
end
hx = xlabel({xlab},"fontsize",12);
ylabel(ylab,"fontsize", 12);
grid on;
title (tit,"fontsize",12);
legend (labels);
