ex10ch2.m

these codes are for solving OED with matlab

function ex10ch2
N = 10;
h = 1/N;
tout = zeros(1,N);
yout = zeros(2,N);
t0 = 0;
y0 = [0; 1];
M = input('Input method: 1 = Heun, 2 = AM2si, 3 = AM2ni ');
tn = t0;
yn = y0;
for n = 1:N
    switch M
    case 1
        [tn,yn] =  Heun(tn,yn,h,@odes);
    case 2
        [tn,yn] = AM2si(tn,yn,h,@odes);
    case 3
        [tn,yn] = AM2ni(tn,yn,h,@odes,@dfdy);
    otherwise
        error('Method not recognized.')
    end
    tout(n) = tn;
    yout(:,n) = yn;
end
tout = [t0 tout];
yout = [y0 yout];
plot(tout,yout,'*',tout,[sin(tout); cos(tout)]);
switch M
case 1
    title('Heun''s method')
case 2
    title('Trapezoidal rule with simple iteration.')
case 3
    title('Trapezoidal rule with Newton iteration.')
end
%======================================================
function dydt = odes(t,y)
dydt = [y(2); -y(1)];

function J = dfdy(t,y)
J = [0 1; -1 0];

function [tnp1,ynp1] = Heun(tn,yn,h,f)
f1 = feval(f,tn,yn);
yn1 = yn + h*f1;
tnp1 = tn + h;
f2 = feval(f,tnp1,yn1);
ynp1 = yn + 0.5*h*(f1 + f2);

function [tnp1,ynp1] = AM2si(tn,yn,h,f)
f1 = feval(f,tn,yn);
p = yn + h*f1;
tnp1 = tn + h;
for i = 1:10
    f2 = feval(f,tnp1,p);
    c = yn + 0.5*h*(f1 + f2);
    if norm(c - p) < 1e-3*norm(c)
        ynp1 = c;
        return
    end
    p = c;
end
error('Simple iteration failed to converge.')

function [tnp1,ynp1] = AM2ni(tn,yn,h,f,dfdy)
f1 = feval(f,tn,yn);
p = yn + h*f1;
tnp1 = tn + h;
J = feval(dfdy,tn,yn);
[L,U] = lu(eye(size(J)) - 0.5*h*J);
for i = 1:10
    f2 = feval(f,tnp1,p);
    residual = (yn + 0.5*h*(f1 + f2)) - p;
    correction = U \ ( L \residual);
    c = p + correction;
    if norm(correction) < 1e-3*norm(c)
        ynp1 = c;
        return
    end
    p = c;
end
error('Newton iteration failed to converge.')