M=2;
T=1;
t = linspace(0,T,M+2);
N = length(t);

U=zeros(M+2,1);
U(1)=1;
for i=1:M+1
   U(i+1)=((t(i+1)^3-t(i)^3)/3+U(i)*(1-2*(t(i+1)-t(i)))/(1+2*(t(i+1)-t(i))));
end
figure
hold on
plot(t,U)
tmore=linspace(0,T,100);
u=(8*tmore.^2-4.*tmore+31*exp(-4.*tmore)+1)/32;

plot(tmore,u,'c')