Why does my plot not show any lines?

%%Variables and values
rhoi=.00139; %m
rhoa=.00199; %m
ri=.0016; %mm
c=11.2; %kPa
c1=0.0499; 
c2=1.0672;
c3=0.4775;
c4=0.0042;
c5=0.0903;
c6=0.0585;
lam=1.8;
Pi=120/7.50061683; %mmHg to kPa
Pa=0; %kPa
rho=rhoi:0.01:rhoa;
e=2.71828;

%%Equations
for r=sqrt((1/lam).*(rho.^2-rhoi.^2)+ri.^2);
F=[rho/r.*lam 0 0;0 r./rho 0;0 0 lam];
C=F.'*F;
I=[1 0 0;0 1 0;0 0 1];
E=(1/2)*(C-I);
Err=E(1,1);
Ett=E(2,2);
Ezz=E(3,3);
Q=(c1*Err^2)+(c2*Ett^2)+(c3*Ezz^2)+(2*c4*Err*Ett)+(2*c5*Ett*Ezz)+(2*c6*Err*Ezz);
trrex=((rho/(r*lam)).^2)*c*e*Q*(2*c1*Err+2*c4*Ett+2*c6*Ezz); %trr extra
tthetaex=((r/rho).^2)*c*e*Q*(2*c2*Ett+2*c4*Err+2*c5*Ezz); %tthetatheta extra
fh=@(r) (1./r)*(tthetaex-trrex);
int=integral(fh,ri,r); %integral to input into lagrange
lagrange=Pi+trrex-int;
dWdErr=c*exp(Q)*((2*c1*Err)+(2*c4*Ett)+(2*c6*Ezz)); %partial of Wfung wrt Erhorho
dWdEtt=c*exp(Q)*((2*c2*Ett)+(2*c4*Err)+(2*c5*Ezz)); %partial of Wfung wrt Ethetatheta
dWdEzz=c*exp(Q)*((2*c3*Ezz)+(2*c5*Ett)+(2*c6*Err)); %partial of Wfung wrt Ezz
trr=-lagrange+((rho/(r*lam))^2)*dWdErr; %cauchy in rr
tthth=-lagrange+((r/rho)^2)*dWdEtt; %cauchy in thetatheta
tzz=-lagrange+(lam^2)*dWdEzz; %cauchy in zz

plot(r, trr,'-')
plot(r, tthth,'-')
plot(r, tzz,'-')

end
%Here is the integral part
function Q = integral(fun,a,b,varargin)
narginchk(3,inf);
if ~isa(fun,'function_handle')
    error(message('MATLAB:integral:funArgNotHandle'));
end
if ~(isscalar(a) && isfloat(a) && isscalar(b) && isfloat(b))
    error(message('MATLAB:integral:invalidEndpoint'));
end
opstruct = integralParseArgs(varargin{:});
Q = integralCalc(fun,a,b,opstruct);

[Star Strider]

Look at your variables.

Note that rho is a single value, not a vector.

I didn’t follow it beyond that, because I have no idea what you want to do.