%Lecture22c
%Integration by Crank-Nicholson Method

clear all

%ODE
%dy/dt=lambda*y
%y(0)=yo

lambda=-1;

%Initial condition
yo=5;

%Integration step
h=input('Enter Integration Step: ');

%Initial time
to=0;

%Final time
tf=input('Enter Final Time: ');

%Integration points
N=(tf-to)/h

i=1;
y(1)=5;
t_vec(1)=0;

for i=1:N
   y(i+1)=y(i)*(2+h*lambda)/(2-h*lambda);
   tvec(i+1)=i*h;
end

%Analytical Solution
tveca=[to:0.001:tf];
ya=yo*exp(lambda*tveca);

figure(1)
plot(tveca,ya,'r-',tvec,y,'b:')
title(['Time Step h=',int2str(floor(h)),'.',int2str((h-floor(h))*10),int2str((h*10-floor(h*10))*10)])
xlabel('t')
legend('analytical','numerical')



