lambda = -2000;
tstar = 1.5;
ts = 40;
k = tstar/ts;
y0 = 0;
yei(:,1) = y0;
ytrap(:,1) = y0;
t = 0;
for n = 1:ts
  yei(:,n+1) = (yei(:,n)-k*lambda*cos(t+k))/(1-k*lambda);
  ytrap(:,n+1) = (ytrap(:,n)+k/2*lambda*(ytrap(:,n)-cos(t))-...
		  k/2*lambda*cos(t+k))/(1-k/2*lambda);
  t = t+k;
end
plot(linspace(0,tstar,ts+1),ytrap,'-*',linspace(0,tstar,ts+1),yei,'-o')
legend('Trapezi','Eulero implicito')
xlabel('t')
ylabel('y(t)')