Double or higher order integration by numerical method

%% For example I am using trapazoidal rule for integration as below

%%

function s=traorl(fun,a,b,M)

h=(b-a)/M

s=0

for k=1:(M-1)

x=a+k*h

s=s+feval(fun,x)

end

s=h*(feval(fun,a)+feval(fun,b))/2+h*s

end

%% I using fun=@(x) a*sin(y)*exp(-b*x+c), and integration with respect to x is done as below

syms a b c x y

fun=@(x) a*sin(y)*exp(-b*x+c)

A=traorl(fun,0,1,5)

%% Integration with respect to x is ok. I wants to integrate again with respect to y by the same function "traorl", But How can I do it as "traorl" function integrate with %%respect to x only. If there is another suggestion or solution do you have then please give me.

%% I wants to do integration only by these numerical method, not by direct integral command.

Answer 1

[Torsten]

Divide the y-interval of integration into subintervals

ystart = y1 < y2 < y3 < ... < yn = yend

and call your function "traorl" in a loop for all these values for y.

Let the results be Iy1,...,Iyn.

Finally call "traorl" with these Iy1,...,Iyn to sum them up in y-direction.

Best wishes

Torsten.

[SHARAD KUMAR UPADHYAY]

Not useful this, please give a short program as .m file