梯形公式matlab程序

function yp=funst(x,y)                      %先建立函数文件
yp=sec(x)+y*tan(x);     

function [t,y]=tixing(f,a,b,y0,N)           %定义梯形公式的函数
h=(b-a)/N;                                  %等距剖分
t=zeros(N+1,1);                             %离散点列向量
y=zeros(N+1,1);                             %数值解列向量
y(1)=y0;                                    %把初值y0附到数值解第一个节点位置
for i=1:N+1                                 %用循环给离散点赋值
 t(i)=a+(i-1)*h;
end
for j=2:N+1                                 %数值解第二个节点开始用迭代求得
       k=1;yc=[];                           %yc是校正时的解数组，k从1开始，迭代一次+1        
       yc(k)=y(j-1)+h*f(t(j-1),y(j-1));     %先用向前Euler公式进行预测
       while k>=1
            k=k+1;
            yc(k)=y(j-1)+h/2*(f(t(j-1),y(j-1))+f(t(j),yc(k-1))); %再用梯形公式校正
            if max(abs(yc(k)-yc(k-1)))<1e-6