%This code works for Problem B in Project 1
%define index vector and constant a
nmax=30; n=(0:nmax)'; a=10; %you may test different inputs
%compute the bounds by hands, vectorize computation
y_u=1./(a*(n+1));
y_l=1./((a+1)*(n+1));

%forwarding recurrence iterations
y_f=zeros(nmax+1,1); %zeros vector for store y
y_f(1)=log((1+a)/a); %initial value
tic
for k=1:nmax
y_f(k+1)=1/k-a*y_f(k); %forward formula
end
t_forward=toc;

%backward recurrence iterations
y_b=zeros(nmax+1,1);
y_b(end)=1/((1+a)*(nmax+1)); %lower bound as initial
tic
for k=nmax:-1:1
y_b(k)=(1/k-y_b(k+1))/a; %backward formula
end
t_backward=toc;

%print out the table using for loop and printf, you may also use table()
fprintf('n\t Lower\t\t Forward\t Upper\t\t Backward\t Quad\n');
y_q=zeros(nmax+1,1);
tic
for k=1:nmax+1
yfun=@(x) x.^(k-1)./(x+a); %define anonymous function
y_q(k)=integral(yfun,0,1); %may also use quad
fprintf('%d\t%f\t %g\t %f\t %f\t %f\n',k-1,y_l(k),y_f(k),y_u(k),y_b(k),y_q(k));
end
t_integral=toc;

fprintf('CPU times: forward=%1.4e s, backward=%1.4e s, integral=%1.4e s\n',t_forward,t_backward,t_integra);
%plots all computed vectors except y_f (too large to fit into same window)
subplot(1,2,1); plot(n,y_l,'--',n, y_u,'--',n, y_b,'-rx', n,y_q,'-.b')
legend('Lower','Upper','Backward','Integral'); set(gca(),'FontSize',16); %set up larger font size
xlabel('n-iterations'); ylabel('y_n');
subplot(1,2,2); plot(n,y_f,'-r.')
legend('Forward'); set(gca(),'FontSize',16); %set up larger font size
xlabel('n-iterations'); ylabel('y_n');