Code gets stuck when generating third set of plots in a for loop

I am writing code for a project which currently involves generating 3 sets of plots (6 plots total) with 1 set per Pr number which there are 3 of (0.1, 1.0, 10). As of now, the current code is capable of generating sets for the first two Pr numbers, but gets stuck indefinitely while trying to generate the plot set for Pr number 10. I am unsure of why this is happening, and any input would be appreciated.

clc
clear
format long 

global Pr 
ddf = 1; 
dtheta = -1;
span = [0 10];

% Increments 
da = 1e-4; 
db = 1e-4;
% Boundary value error tolerance
tol = 1e-6;
% Maximum number of iterations
itmax = 1000;
% Prandtl numbers
Prm = [0.1, 1.0, 10];
for k = 1:2
    % Primary shot
    Pr = Prm(k);
    for i = 1:itmax
        init = [0, 0, ddf, 1, dtheta];
        % obtain boundary values of f' and theta at eta = eta max
        [e, z] = ode45(@odefun, span, init); 
        % Store boundary values at eta max
        df_k = z(end,2);
        theta_k = z(end, 4);

        if df_k < tol && theta_k < tol
            break;
        else
            for j = 1:2
                if j == 1
                    init2 = [0, 0, ddf+da, 1, dtheta];
                    [e, z1] = ode45(@odefun, span, init2);
                    df_ka = z1(end,2);
                    theta_ka = z1(end,4);

                else
                    init3 = [0, 0, ddf, 1, dtheta+db];
                    [e, z2] = ode45(@odefun, span, init3);
                    df_kb = z2(end,2);
                    theta_kb = z2(end,4);
                end
            end
            fa = (df_ka-df_k)/da;
            fb = (df_kb-df_k)/db;
            thetaa = (theta_ka-theta_k)/da;
            thetab = (theta_kb-theta_k)/db;
            A = [fa fb; thetaa thetab];
            b = [-df_k; -theta_k];
            v = A\b;
        end
        ddf = ddf + v(1);
        dtheta = dtheta + v(2);
    end
    figure
    hold on
    plot(e, z(:,2));
    plot(e, z(:,4));
    yline(0)
    xlabel('{\eta}')
    ylabel('df/d{\eta} & {\theta}')
    legend('df/d{\eta}', '{\theta}')

    figure
    hold on
    plot(e, z(:,3))
    plot(e, -z(:,5))
    yline(0)
    xlabel("{\eta}")
    ylabel("d^2f/d{\eta}^2 & -d{\theta}/d{\eta}")
    legend("d^2f/d{\eta}^2", "-d{\theta}/d{\eta}")
end

%% ODE Functions
function dfde = odefun(~, z)
global Pr
dfde = [ z(2)
    z(3)
    2*(z(2)^2)-3*z(1)*z(3)-z(4)
    z(5)
    -3*Pr*z(1)*z(5) ];
end

Answer 1

[VBBV]

clc
clear
format long 

% global Pr 
ddf = 1; 
dtheta = -1;
span = [0 10];

% Increments 
da = 1e-4; 
db = 1e-4;
% Boundary value error tolerance
tol = 1e-6;
% Maximum number of iterations
itmax = 10000;
% Prandtl numbers
Prm = [0.1, 1.0, 10];
for k = 1:3
    % Primary shot
    Pr = Prm(k);
    for i = 1:itmax
        init = [0, 0, ddf, 1, dtheta];
        % obtain boundary values of f' and theta at eta = eta max
        [e, z] = ode45(@odefun, span, init); 
        % Store boundary values at eta max
        df_k = z(end,2);
        theta_k = z(end, 4);

        if df_k < tol && theta_k < tol
            break;
        else
            for j = 1:2
                if j == 1
                    init2 = [0, 0, ddf+da, 1, dtheta];
                    [e, z1] = ode45(@odefun, span, init2);
                    df_ka = z1(end,2);
                    theta_ka = z1(end,4);

                else
                    init3 = [0, 0, ddf, 1, dtheta+db];
                    [e, z2] = ode45(@odefun, span, init3);
                    df_kb = z2(end,2);
                    theta_kb = z2(end,4);
                end
            end
            fa = (df_ka-df_k)/da;
            fb = (df_kb-df_k)/db;
            thetaa = (theta_ka-theta_k)/da;
            thetab = (theta_kb-theta_k)/db;
            A = [fa fb; thetaa thetab];
            b = [-df_k; -theta_k];
            v = A\b;
        end
        ddf = ddf + v(1);
        dtheta = dtheta + v(2);
    end
    figure
    hold on
    plot(e, z(:,2));
    plot(e, z(:,4));
    yline(0)
    xlabel('{\eta}')
    ylabel('df/d{\eta} & {\theta}')
    legend('df/d{\eta}', '{\theta}')

    figure
    hold on
    plot(e, z(:,3))
    plot(e, -z(:,5))
    yline(0)
    xlabel("{\eta}")
    ylabel("d^2f/d{\eta}^2 & -d{\theta}/d{\eta}")
    legend("d^2f/d{\eta}^2", "-d{\theta}/d{\eta}")
end

%% ODE Functions
function dfde = odefun(Pr, z)
% global Pr
dfde = [ z(2)
    z(3)
    2*(z(2)^2)-3*z(1)*z(3)-z(4)
    z(5)
    -3*Pr*z(1)*z(5) ];
end