I solved the problem in another approach. I had the points on the function, so I implemented the trapezoidal rule, though I'd like to know what I should have fixed in the code (more cleaner now).
function [C_l]=lift_coefficient(a,b) % the function calculates the lift coeeficient of a wing's cross-sections %inputs - vectors a and b are defining the cross-sectional parameters
if nargin<2 a=[0.3 0.5 0.7]; b=[-0.7 -0.8 -0.9]; end;
C_l=quad(@camber,-0.9999,0.9999) %numerical integration of f
function f=camber(x) y_c=a./2*(1-x.^2)./(2+b.*(x+1)); %camber-line function dyc_dx=diff(y_c,x)% first derivative of y_c f=-2*sqrt((1+x)./(1-x))*dyc_dx % integrand end
end