% Exact and numerical solutions to Problem 2.30 of the text
%
%    2
%   d y(t)     dy(t)
%   ------ + 2 -----  + y(t) = x(t),  where x(t) = 0
%      2
%    dt         dt
%
%  y(0) = 2, y'(0) = -1 
%  


clf

t = 0:0.001:10;
y = 2*exp(-t) + t.*exp(-t);
plot(t,y);
xlabel('t (seconds)');
ylabel('y(t)');
grid
title('Exact and Numerical Solutions to d^2y(t)/dt^2 + 2 dy(t)/dt + y(t) = 0');
hold on

%
% y[n+2] - 2 y[n+1] + y[n]      y[n+1] - y[n]
% ------------------------  + 2 ------------- + y[n] = x[n]
%            2
%           T                         T
%
%                                2         2
% y[n+2] + (2T-2)y[n+1] + (1-2T+T )y[n] = T x[n]
%
% y(0) = 2 => y[0] = 2
%
%          y[n+1] - y[n]
% y'(t) =  -------------           
%               T
%
%          y[1] - y[0]
% y'(0) =  ----------- = -1    => y[1] = y[0] - T = 2 - T
%               T

T = 0.4;
a = [2*T-2 1-2*T+T^2]; % a1, a2
b = [0 0 T^2];         % b0, b1, b2
x0 = [0 0];            % x[0] x[1]
y0 = [2 2-T];          % y[0] y[1]
n = 0:10/T;            % Go out to 10 seconds
x = zeros(1,length(n));
y1 = recur(a,b,n,x,x0,y0);
plot(n*T,y1,'o');

T = 0.1;
a = [2*T-2 1-2*T+T^2]; % a1, a2
b = [0 0];             % b0, b1, b2
x0 = 0;                % x[0] x[1]
y0 = [2 2+T*(-1)];     % y[0] y[1]
n = 0:10/T;            % Go out to 10 seconds
x = zeros(1,length(n));
y2 = recur(a,b,n,x,x0,y0);
plot(n*T,y2,'sr');
legend('Exact','T = 0.4','T = 0.1');
print -dpsc2 p2_30.ps
