% Igor Yanovsky
% Teaching Assistant
% Math 151B
% 04/15/07
%
% Homework 2, Problem 4(c)
%
%#####################################################################
%
% This program solves the ordinary differential equation
%
% dy/dt = 1+y/t, 1 <= t <= 2
% y(1) = 2
%
% using RUNGE-KUTTA MIDPOINT 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

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

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

t = t0;
w(1) = y0; % initial condition
yexact(1) = y0;
error(1) = 0; % error is 0 initially

for(i=1:N)
w(i+1) = w(i) + h*f( t + (h/2), w(i)+(h/2)*f(t,w(i))); % advances the ODE using RK MIDPOINT method
t = t + h;
yexact(i+1) = findExact(t);
end

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

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(N+1)); fprintf(1,',\n');
fprintf(1,'error = %5.10e', error ); fprintf(1, '\n');

plot(yexact,'r','LineWidth',2)
hold on
plot(w,'b--','LineWidth',1.5)
legend('exact','Midpoint')

clear all