function [t,w]=imptrap(func,diffunc,a,b,n,alpha) h=(b-a)/n; t=linspace(a,b,n+1)'; w=zeros(n+1,1); w(1)=alpha; for i=1:n old=w(i); err=1; j=1; while err>1e-12 j=j+1; new=old-(old-w(i)-h/2*(func(t(i),w(i))+func(t(i+1),old))/... (1-h/2*diffunc(t(i+1),old))) err=abs(new-old)/old; old=new; if j>100, error('max iterations exceeded'); end end w(i+1)=old; end end