% Igor Yanovsky
% Teaching Assistant
% Math 151B
% 04/15/07
%
% Homework 2, Problem 4(b)
%
%#####################################################################
%
% This program solves the ordinary differential equation
% 
%  dy/dt = 1+y/t,    1 <= t <= 2
%  y(1) = 2
% 
% using Euler's method. The function defining the ODE is defined in a
% separate Matlab function.
%
% This program also finds the error of the approximation by comparing the
% approximate solution and the exact solution  y(t) = t log(t) + 2t.
% 
%#####################################################################

clear all

format long

% ODE Specification
%

t0 = 1.0;
tF = 2.0;
y0 = 2.0;

h = 1.164008846467233e-004   % timestep
N = (tF-t0)/h                % number of timesteps

tn = t0;
wn(1)     = y0;              % initial condition
yexact(1) = y0;

for(i=1:N)
    wn(i+1) = wn(i) + h*f(wn(i),tn); % advances the ODE one step using Euler's method
    tn = tn + h;
    yexact(i+1) = tn * log(tn) + 2*tn;
end

tn;

error = abs(wn(N+1) - yexact(N+1));

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

plot(wn,'b')

plot(yexact,'r')
   
clear all