The problem is that you have many values of lambda. To deal with this use a loop to do the integration for each value of lambda (this requires f1 to have an extra input parameter):
lambda=1:0.1:2;
gamma1=0.4;
zeta=0.3;
a=0;%lower limit
b=1;%lower limit
n=10;%number of sub-intervals
f1=@(x,p)((1-x).^((gamma1./zeta)-1)).*(exp(-(lambda(p).*x)./zeta));
h=(b-a)/n;
for p = 1:numel(lambda)
for k=1:n
x(k)=a+k*h;
y(k)=f1(x(k),p);
end
so=0;se=0;
for k=1:n-1
if rem(k,2)==1
so=so+y(k);%sum of odd terms
else
se=se+y(k); %sum of even terms
end
end
SOL(p)=h/3*(f1(a,p)+f1(b,p)+4*so+2*se);
end
