echo off;
% NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
% To accompany the text:
% NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
% Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
% This free software is complements of the author.

% Algorithm 9.8 (The Hamming Method).
% Section 9.6, Predictor-Corrector Method, Page 473
echo on; clc; format long; hold off; clear

% This program implements Hamming's method

% for solving the initial value problem

% y' = f(t,y) with y(a) = y0

% Define and store the function f(t,y) in the M-file f.m
% function z = f(t,y)
% z = 30 - 5*y;

delete f.m
diary f.m; disp('function z = f(t,y)');...
disp('z = 30 - 5*y;');...
diary off;

f(0,0); % Remark. f.m rk4.m hamming.m are usd for Algorithm 9.8
pause % Press any key to continue.

clc;
% Place the endpoints of [a,b] in a and b.

% Place the initial value y(a) in ya.

% Place the number of subintervals in n.

a = 0;
b = 10;
ya = 1;
n = 70;

pause % Press any key to continue.

clc;
h = (b - a)/n;
T = zeros(1,n+1);
Y = zeros(1,n+1);

% Generate three starting values using the Runge-Kutta method.
[Ts,Ys] = rk4('f',a,a+3*h,ya,6);
T(1:4) = Ts(1:2:7);
Y(1:4) = Ys(1:2:7);

% Proceed with Hamming's method.
[T,Y] = hamming('f',T,Y);
points = [T;Y];

pause % Press any key to see the list of solution points.

clc;, clg;
c = 0;
d = 7;
axis([a b c d]);...
plot(T,Y,'-r');...
hold on;...
plot([a b],[0 0],[0 0],[c d]);...
xlabel('t');...
ylabel('y');...
Mx1 = 'Hamming`s solution to y` = f(t,y).';...
title(Mx1);...
grid;...
shg; pause % Press any key to continue.

Mx2 = ' t(k) y(k)';
Mx3 = 'Hamming`s method is stable';
Mx4 = ['for n = ',num2str(n),' and h = '];
Mx5 = [Mx4,num2str(h),' because'];
Mx6 = ' h < 0.69/|f (t,y)|';
Mx7 = ' y';
clc,echo off,diary output,...
disp(''),disp(Mx1),disp(''),disp(Mx2),disp(points'),...
disp(Mx3),disp(Mx5),disp(''),disp(Mx6),disp(Mx7),...
diary off, echo on

pause % Press any key to continue.

clc;
% Place the endpoints of [a,b] in a and b.

% Place the initial value y(a) in ya.

% Place the number of subintervals in n.

a = 0;
b = 10;
ya = 1;
n = 50;

pause % Press any key to continue.

clc;
h = (b - a)/n;
T = zeros(1,n+1);
Y = zeros(1,n+1);

% Generate three starting values using the Runge-Kutta method.
[Ts,Ys] = rk4('f',a,a+3*h,ya,6);
T(1:4) = Ts(1:2:7);
Y(1:4) = Ys(1:2:7);

% Proceed with Hamming's method.
[T,Y] = hamming('f',T,Y);
points = [T;Y];

pause % Press any key to see the list of solution points.

clc;
plot(T,Y,'-g');...
xlabel('t');...
ylabel('y');...
title(Mx1);...
grid;...
axis;...
hold off;...
shg; pause % Press any key to continue.

Mx3 = 'Hamming`s method is unstable';
Mx4 = ['for n = ',num2str(n),' and h = '];
Mx5 = [Mx4,num2str(h),' because'];
Mx6 = ' 0.69/|f (t,y)| < h';
Mx7 = ' y';
clc,echo off,diary output,...
disp(''),disp(Mx1),disp(''),disp(Mx2),disp(points'),...
disp(Mx3),disp(Mx5),disp(''),disp(Mx6),disp(Mx7),...
diary off, echo on