来自「此文挡是用Matlab开发的几种计算数值微积分的方法」· 代码 · 共 26 行

part1
function s=trapr1(f,a,b,M)
h=(b-a)/M;
s=0;
for k=1:(M-1)
    x=a+h*k;
    s=s+feval(f,x);
end
s=h*(feval(f,a)+feval(f,b))/2+h*s;
***************************
part2
function y1=f(x)
y1=exp(x);
***************************
test
>> s=trapr1('f',0,1,200)

s =

    1.7183

>> vpa(s)
 
ans =
 
1.71829