function [t,y]=eulertp(f,tspan,startval,step,K)
% the only difference between this and feuler is the last input
% K is the max number of iterates for the while loop of the iterative
% method based on trapezoidal rule.

t=[tspan(1):step:tspan(2)];
steps=length(t);
y=zeros(1,steps);
y(1)=startval;

for i=2:steps
    y1=y(i-1)+step*feval(f,t(i-1),y(i-1));
    
    % notice this is not yet y(i), y1 is only an intermediate geuss at y(i) 
    % which we'll improve by an iterative method based on trapezoidal rule
    it=0;
    diff=1;
    while (diff>.05)&(it<K)
        y2=y(i-1) + step * (feval(f,t(i-1),y(i-1)) + feval(f,t(i),y1))/2;
        diff=abs(y2-y1);
        y1=y2;
        it=it+1;
    end
    
    % now set y(i) to the final value from the iterative method
    y(i)=y1;
end