function [ t,y ] = pece2(f,tspan,y_0,N_s,yes) 

% solve the ODE dy/dt = f(t,y) by 2nd-order PECE method in N_s steps,
% with midpoint rule as predictor and trapezoid rule as corrector 

t_0=tspan(1);
t_f=tspan(2);
D=length(y_0);

dt = (t_f - t_0)/N_s;

t = t_0:dt:t_f;
N=length(t);

yj=y_0';
y(1,:) = y_0;
j = 1; 

% first step by Euler method 
yd(1,:)=feval(f,t(1),yj)';
neval=1;
yj = yj + dt*yd(1,:)';
y(2,:) = yj';
j=2;

while j < N 
yj1=y(j-1,:)';
yj0=y(j,:)';
fj0=feval(f,t(j),yj0);
neval=neval+1;
yj = yj1 + 2*dt*fj0;
yj = yj0 + dt*(fj0+feval(f,t(j+1),yj))/2;
neval=neval+1;
y(j+1,:) = yj';
j = j + 1;
end

if yes==1

for k=1:D,

figure
z=y(:,k);
plot(t,z) 
xlabel('time t')
ylabel(sprintf('y_%d', k))


end

end 


t=t'; 

% neval

return