ex14ch2.m

these codes are for solving OED with matlab

function ex14ch2
% plot the boundary of the absolute stability region for
% five linear multistep methods
for method = 1:5
    ntheta = 180;
    theta1 = 0;
    theta2 = pi;
    theta3 = (theta2 - theta1) / ntheta;
    i = sqrt(-1);
    index = 0;
    u = zeros(ntheta,1);
    v = zeros(ntheta,1);
    for itheta = 1:ntheta+1
        theta = theta1 + (itheta-1)*theta3;
        r = exp(i * theta);
        % Forward Euler AB1:
        if method == 1
            rho = -r + 1;
            sigma = 1;
        end
        % Backward Euler BDF1:
        if method == 2
            rho = -r + 1;
            sigma = r;
        end
        if method == 3
            rho = -r^2 + (4/3) * r - 1/3;
            sigma = (2/3) * r^2;
        end
        % Adams-Moulton AM3:
        if method == 4
            rho = -r^2 + r;
            sigma = (5/12)*r^2 + (8/12)*r - (1/12);
        end
        % Backward Differentiation BDF3:
        if method == 5
            rho = -r^3 + (18/11)*r^2 - (9/11)*r + (2/11);
            sigma = (6/11)*r^3;
        end
        if sigma ~= 0
            index = index + 1;
            z = - rho / sigma;
            u(index) = real(z);
            v(index) = imag(z);
        end
    end
    figure
    plot(u,v,u,-v)
    axis on
    xlabel('Real (h\lambda)')
    ylabel('Imag (h\lambda)')
    title(['Absolute Stability Region Boundary for Method ',int2str(method)])
end