Non-divergence as boundary condition for ODE

I have a system of two second-order ODEs.

- The ODE for is:

where:

- The ODE for is:

where:

where the coefficients and depend on the the first derivative and where the parameters are given.

% Parameters
gamma = 4;
theta2 = 0.0025;
theta1 = - 0.0150;
sigmaD = 0.0240;
r = 0.0041;
delta = 1;
p12 = 0.1000;
p21 = 0.0167;
pi2 = p12 / (p12 + p21);
Gamma_pi = (theta2 - theta1) / (r * (r + p12 + p21));

I want to solve this system imposing as sole boundary condition on the ODEs of and the fact that they must not diverge to either or over the interval . I start by writing the system of second order ODEs as a sytem of first order ODEs in an ode_Sf.m file:

function dy = ode_Sf(pi, y, gamma, theta2, theta1, sigmaD, r, delta, p12, p21, pi2, Gamma_pi)
    
    h = (theta2 - theta1) / sigmaD * pi * (1 - pi);

    Q3 = h^2 / 2;
    Q1 = gamma * sigmaD * h + r * gamma * Gamma_pi * h^2 - (p12 + p21) * (pi2 - pi);
    Q0 = (gamma * r)^2 / 2 * Gamma_pi^2 * h^2 + gamma^2 * r * Gamma_pi * sigmaD * h + r * log(delta);

    P3 = Q3;
    P1 = gamma * r * Gamma_pi * h^2 + gamma * sigmaD * h - (p12 + p21) * (pi2 - pi) + y(2) * h^2;
    P0 = gamma * r * Gamma_pi^2 * h^2 + 2 * gamma * Gamma_pi * sigmaD * h + y(2) * Gamma_pi * h^2 + y(2) / r * sigmaD * h;

    dy = [y(2)
          y(2)^2 + y(2) * Q1 / Q3 + y(1) * r / Q3 + Q0 / Q3
          y(4)
          y(4) * P1 / P3 + y(3) * r / P3 + P0 / P3];
end

Then I define the model and find the solution over where the interval is smaller than since the ODEs have singular points at and .

%% Model
pi_range = linspace(0.001, 0.995, 994);
y0 = [-0.0001 -300 -105 -30];
model = @(pi, y) ode_Sf(pi, y, gamma, theta2, theta1, sigmaD, r, delta, p12, p21, pi2, Gamma_pi);
[pi, y] = ode45(model, pi_range, y0);

Warning: Failure at t=9.889078e-01. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.776357e-15) at time t.

I want my final solutions to be negative, convex and U-shaped. I provide here an example for . Notice that in my code I set while my final implementation should include different values of gamma.

Unfortunately, my solutions explode to infinity and I don't know how to solve this problem.

% Plot f
figure;
plot(y(:,1));
xlabel('\pi');
ylabel('f(\pi)');
grid on;

% Plot S
figure;
plot(y(:,3));
xlabel('\pi');
ylabel('S(\pi)');
grid on;

Could you advise how to impose the boundary condition that the ODEs do not explode?