29.03.2016

lab_6.m

clear;
% Initializing ode45 options 
options = odeset('AbsTol', 1.e-9, 'RelTol', 1.e-10);

f=@(x,u)[u(2);-u(3)*u(1);0];

% values to test:
% [L1, L2] = [0.5, 0.6],
% [L1, L2] = [2.2, 2.3] 
% [L1, L2] = [, ],

values = [0.5, 0.6; 2.2, 2.3; 8.0, 10.0];

for value = values',

    L1 = value(1); 
    sol = ode45(f, [0, pi], [0, 1, L1], options);
    U1 = sol.y(1, end);

    L2 = value(2); 
    sol = ode45(f, [0, pi], [0, 1, L2], options);
    U2 = sol.y(1, end);

    err = 1;
    eps = 1.e-8;
    k = 1;
    k_max = 100;
    
    while err > eps && k < k_max
      k = k + 1;
    
        L = L2 - U2 * (L2 - L1)/(U2 - U1);
        sol = ode45(f, [0, pi], [0, 1, L], options);
        U = sol.y(1, end);
    
        U1 = U2;
        U2 = U;
        L1 = L2;
        L2 = L;
    
        err = abs(U);
    end

    plot(sol.x, sol.y(1,:))
    legend(['U', num2str(L)],)
    hold on
end