clear all;

Nt=30;

T=linspace(0,1,Nt);

y0=1;

Y1=zeros(1,Nt);
Y2=Y1;
Y3=Y1;



Y1(1)=y0;Y2(1)=y0;Y3(1)=y0;
delta=T(2)-T(1);


P(1) = exp( coeff(T(1)) )*delta ;
for it = 2:Nt
    P(it) = exp( (coeff(T(it))+coeff(T(it-1) ))*delta/2 );
Y1(it) = P(it)*Y1(it-1);

I = T(it)^3/3+T(it) - T(it-1)^3/3 - T(it-1);
Y2(it) = exp(I)*Y2(it-1);

    
Y3(it) = exp(T(it)^3/3+T(it))*y0;

end

figure(1)
subplot(2,1,1)
plot(T,Y1,'*b',T,Y2,'or',T,Y3,'k')
subplot(2,1,2)
plot(T(2:end),P(2:end))