function S=trapezia(f,a,b,f0,h)
%f為函數f(x,y);
%a,b為取值范圍的左右端點；
%f0為函數初值y(0)，h為步長,n為弟迭代步數；
%feval函數實現輸入變量按照指定函數計算的函數值.
n=(b-a)/h;
T=zeros(1,n+1);
Y=zeros(1,n+1);
T=a:h:b;
Y(1)=f0;
for j=1:n
    z0=Y(j)+h*feval('f',T(j),Y(j));
    Y(j+1)=Y(j)+h*(feval('f',T(j),Y(j))+feval('f',T(j+1),z0))/2;
end
t=[0:.1:5];
y=1./exp(10*t);
figure(1)
plot(T,Y,'*',T,Y,'r-',t,y,'k-');
hold on
T=[T' Y']