This
bl = @(ln1, ln2) bl+pro;
is syntactically wrong. The operation + is not defined on function handles. You need to include the arguments. The statement should read
bl = @(ln1,ln2) bl(ln1,ln2) + pro(ln1,ln2);
You'll need to make similar changes in the definition of cpd and pro. Try some simple examples before jumping to the real thing to make sure you understand how the language works first. Something like:
>> f = @(x)ones(size(x));
>> f = @(x)f(x).*x + 1; % now f(x) = x + 1
>> f = @(x)f(x).*x + 1; % now f(x) = x^2 + x + 1
>> f = @(x)f(x).*x + 1; % now f(x) = x^3 + x^2 + x + 1
>> f(1)
ans =
4
The next problem you will encounter here is that you are using matrix multiplication * and matrix division (linear system solves) / when what you really mean is element-wise multiplication .* and element-wise division ./. This will become important when you call integral2 because integral2 requires that your function be able to accept array input and operate element-wise. You are doing this at the last step of defining bl, but it is just as important in the expressions in the loop as when finishing off the definition. Just FYI, if you have a function that can't take an array, you can always integrate
@(x,y)arrayfun(f,x,y)
instead of integrating just f. ARRAYFUN takes care of calling f element-wise with only scalar inputs. I don't think you'll need that here, however, as it looks like everything can be done element-wise.
