Paste: newtons method for stiff DEs

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