function HW12014Ex1(N)
%
% Solves the ivp y_{n+2}+  y_{n} = n^2
%  y_0=y1=1
% 
%
fprintf('\n')
disp('______________________________________________')
disp('i         y_n         yclosed            ')
disp('______________________________________________')
%
ye=zeros(1,N);
yr=zeros(1,N);
yex=zeros(1,N);
ye(1,1:N)=(1/2)*(i).^(0:N-1).*(1+(-1).^(0:N-1))+(3/4)*(i).^(1:N).*(-1+(-1).^(0:N-1))+(1/2)*(0:N-1).^2-(0:N-1);
yex(1,1:N)=cos((0:N-1)*pi/2)+(3/2)*sin((0:N-1)*pi/2)+(1/2)*(0:N-1).^2-(0:N-1);
%
yr(1,1)=1;yr(1,2)=1;
for k=3:N
    y0=yr(1,k-2);y1=yr(1,k-1);
    yr(1,k)=-y0+(k-3).^2;
end
%
for k=1:N
fprintf(' %5.0f  %15.6e %15.6e %15.6e\n',k-1,ye(1,k),yex(1,k),yr(1,k))
end
xx=zeros(1,N);
xx(1,1:N)=(0:N-1);
plot(xx,ye,'r*',xx,yr,'b+')
xlabel('n','Fontsize',14);ylabel('y_n','Fontsize',14);
legend('Recurrence','Closed form','Location','SouthEast')
