% Igor Yanovsky
% Teaching Assistant
% Math 151B
% 04/18/07
%
%#####################################################################
%
% This program solves the ordinary differential equation
% 
%  dy/dt = f(t,y)
%  y(t_0) = y_0
% 
% using IMPLICIT TRAPEZOIDAL method. The function defining the ODE is
% defined in a separate Matlab function.
%
%#####################################################################

clear all

format long

eps = 10e-7;

t0 = 0.0;
tF = 2.0;
y0 = -1.0;

N = 100;     % number of timesteps

h = (tF - t0)/N;

t = t0;
w(1) = y0;
yexact(1) = y0;
error(1) = 0;

for(i=1:N)
    wkm1 = w(i);    %  initial condition for Newton's method
    
    for(k=1:100)    % Newton's iteration
        wk = wkm1 - (wkm1 - w(i) - (h/2.0)*(f(t,w(i))+f(t+h,wkm1) ) / (1-(h/2.0)*fPrime(t+h,wkm1)));
        if( abs(wk-wkm1) < eps )
            break;
        end
        wkm1 = wk;
    end
    
    % advances the ODE one step using Implicit Trapezoidal Method:
    w(i+1) = w(i) + (h/2.0)*( f(t,w(i)) + f(t+h,wk) );
    t = t + h;
end

yexact = findExact( t );
error = w(i+1) - yexact;

fprintf(1,'N = %3d',N);  fprintf(1,', h = %4.7e',h);  fprintf(1,', t = %4.2f',t);  fprintf(1,',\n');
fprintf(1,'w = %5.10e', w(N+1)    );  fprintf(1,',\n');
fprintf(1,'y = %5.10e', yexact);  fprintf(1,',\n');
fprintf(1,'error = %5.10e', error );  fprintf(1, '\n');


clear all